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A  simplified  method  for  predicting  rotor  blade  airloads 
Wang  Shicun  and  Xu  Zhi  (Nanjing  Aeronautical  Institute) 

At  present,  a  simplified  method  to  predict  the  airloads  of  rotor 
blades  is  urgently  needed  for  engineering  applications.  In  this  paper, 
on  the  basis  of  classical  generalized  vortex  theory  of  the  rotor, 
through  the  simplification  process  of  the  contribution  of  the  circula¬ 
tions  to  the  induced  velocities,  the  second  order  distribution  of  the 
induced  velocities  of  the  rotor  is  obtained.  Then,  based  on  the  blade 
element  theory,  a  closed  form  of  equations  for  circulation  is  established. 
By  taking  the  flapping  condition  into  account,  simplified  formulae  for 
predicting  rotor  blade  airloads  are  obtained.  In  particular,  in  the 
expressions  of  flapping  coefficients,  the  variation  of  the  induced 
velocity  is  taken  into  consideration  by  directly  relating  to  blade 
parameters  and  flight  parameters.  Finally,  the  example  given  indicated 
that  these  equations  are  suitable  for  aerodynamic  analysis  and  prelim¬ 
inary  design  of  helicopters. 

SYMBOLS 


(2 — rotor  rotation  velocity 
R — rotor  radius 
» — induced  velocity 
i—v/QR — dimensionless  v 
r — circulation 
r-T/Q/?* — dimensionless  r 

(r,  iji)  or  (p,  0) — polar  coordinate  on  the  plane  of  the  rotor  blade 
r  ■  r/R — dimensionless  r 
p  ■  p/R — dimensionless  p 
PH — density  of  air 
k — number  of  blade 

l - 

This  paper  was  presented  at  the  7th  European  Rotor  Flight  Vehicle  and 
Power  Lift  Flight  Vehicle  Conference  in  September  1981.  Reed.  Nov.  1981. 
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V — forward  velocity  of  the  rotor 
o 

V  ■  V  /ftR — dimensionless  V 
oo  o 

V, — combined  air  flow  velocity,  a  constant  with 

x  respect  to  the  blade  plane 

*  V^/fiR — dimensionless 

aQ — attack  angle  of  rotor  with  respect  to  VQ 
— attack  angle  of  rotor  with  respect  to 
b —  arc  length  of  the  blade 
5  =  b/R — dimensionless  b 

c — lift  coefficient  in  the  cross-section  of  the  blade 

ao — slope  of  two  dimensional  lift  curve 

T1 — pulling  force  of  a  blade 

T  ®  k^ — pulling  force  of  the  rotor 

— installation  angle  of  the  blade  cross  section 

— installation  angle  of  the  cross  section  of  the 
blade  root 

Afl — twist  of  the  blade 

— cosine  term  of  the  periodic  variation  of  blade 
displacement 

— sine  term  of  the  periodic  variation  of  blade 
displacement 

U — relative  flow  velocity  of  the  blade  cross  section 
U  *  U/flR — dimensionless  U 

U — velocity  component  of  U  in  the  plane  of  blades 

A 

U — velocity  component  of  U  in  the  direction 
y  perpendicular  to  the  blades 

cosa;  /flR — advance  ratio 
o  o 

VQ  sinao/8R — flowing-in  ratio 

K — flapping  adjustment  coefficient 
8e — flapping  angle  with  respect  to  the  horizontal  hinge 
0 — flapping  angle  with  respect  to  the  rotation  center 
a0 — coning  angle 

an — coslne  term  of  blade  flapping 
bn — sine  term  of  blade  flapping 
e — outward  displacement  of  the  horizontal  hinge 
e  ■  e/R — dimensionless  e 
m, — mass  of  the  blade 


»i"«»/Pj *R* — dimensionless 

C^«7y-|-p**/?*Q1i?, — coefficient  of  pulling  force 

J  — moment  of  inertia  of  a  blade  around  the  horizon- 
e  tal  hinge 

S — moment  of  mass  of  a  blade  around  the  horizontal 
e  hinge 

(M.) — torque  of  pulling  force  of  a  blade  around  the 
e  horizontal  hinge 

(Mr)  — gravity  torque  of  blade  pulling  force  around  hor- 
e  ijzontal  hinge 

g — gravitational  acceleration  . 

x~’££pe§Pt!2  £?Heioefficient  of  the  root  and  the 

I .  INTRODUCTION 

The  estimation  of  airloads  of  rotor  blades  in  the  flapping  plane 
is  one  of  the  basic  topics  In  the  aerodynamics  and  helicopter  aero¬ 
dynamics.  This  is  because:  the  flight  quality,  pilot  control  quality 
of  the  helicopter  and  the  fatigue  life  of  the  rotor  and  the  destabiliza¬ 
tion  of  the  aerodynamic  elasticity  are  entirely  determined  by  the  air¬ 
loads  of  the  rotor.  It  Is  particularly  important  with  regard  to  its 
strain  load. 

Since  the  1960's,  scholars  In  all  the  countries  have  made  signi¬ 
ficant  progress  In  this  area  through  hard  work.  AGARD  held  a  meeting 
on  "Methods  for  Estimation  of  Helicopter  Rotor  Loads"  in  Italy  in  1973 
[1]  and  demonstrated  the  various  analytical  methods  adopted  by  various 
manufacturers.  However,  Just  as  what  is  described  in  some  comments  and 
papers  published  later  [2,3],  due  to  the  complexity  of  the  working  con¬ 
dition  of  the  rotor,  even  with  a  high  speed  large  scale  numerical  com¬ 
puter,  no  significant  progress  has  been  made  in  recent  years.  Therefore, 
it  is  more  attractive  to  develop  a  simplified  method  to  estimate  rotor 
airloads  for  use  by  design  engineers. 

This  paper  began  with  the  generalized  vortex  theory  of  the  rotor 
to  obtain  the  relations  between  first  two  orders  of  harmonic  wave  of 
the  induced  velocity  and  the  lower  and  same  order  harmonic  of  the  cir¬ 
culations.  Then,  based  on  the  blade  element  theory,  a  closed  form 
equation  for  circulations  are  established.  Finally,  by  taking  flapping 
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conditions  into  account,  the  simplified  equation  for  the  estimation 
of  the  airload  of  the  rotor  is  derived. 

II.  INDUCED  VELOCITIES 

Based  on  the  generalized  vortex  theory  of  rotor  with  fixed  vortex 
[4],  the  induced  velocity  at  any  point  on  the  plane  of  the  blades 
(r,  ty)  is  a  function  of  the  adhered  circulation  r  (p,  9): 

7— 7(r)  (2-1) 

If  the  circulation  r  is  expanded  into  a  Fourier  series: 

Cr«*(p)cosm0+f«#(7)ai»m83  (2-2) 

_  m»l  v  ' 

Then,  the  induced  velocities  v  can  also  be  written  as  a  Fourier  series: 

+  2  C®~(7)co8i4+o„(?)sin,i£i  (2-3) 

Here  the  different  harmonic  wave  components  of  the  induced  velocity 

in  general  are  caused  by  the  various  orders  of  harmonic  components  of 
the  circulations.  In  this  paper,  for  simplification,  only  the  major 
contributing  lower  order  and  same  order  terms  of  the  circulations  are 
considered  with  regard  to  induced  velocities,  i.e. 

5,-55 

«t.-o!.  +  o§I  +  5iJ+53: 
o>/-o|/+®ii+5i;+5|; 


where  the  superscripts  represent  the  order  of  harmonic  of  the  circula¬ 
tions. 

According  to  Wang  Shicun’s  vortex  theory  1  by  limiting  to  second 
order  harmonics,  we  get  (the  downward  induced  velocity  is  positive)  [5]: 


4«Vi 

-b  -H-b  -b  *•  -s> 


♦Jb£ -5-  -bK-b  -b  *•  -£>]+<*-*><-« } 
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where  the  super-geometric  function  is  defined  as: 
where  | z |  <1 

(a)*-  a  (a  +l)(o  +2 )•••(«  +  k  — 1) 

:•  d**-l,  -2,  — 

and  the  symbol  c  is 

e  —  1  —  Isinct,!  _  cos  a, 


cos  a, 


1  +|sma,| 


<1 


(2-5) 


The  above  expressions  for  the  harmonic  components  of  the  induced 
velocities  are  the  key  to  the  solution  in  this  paper.  The  super-geo- 
metric  function  used,  based  on  [8],  can  take  only  the  first  few  terms. 

III.  CIRCULATION  EQUATION 


Based  on  the  blade  element  theory  and  the  famous  Zukovski  equation, 
the  airload  of  the  rotor  blade  can  be  expressed  as  [6]: 
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0 


\  V/.  % 


pi 


(3-D 


(3-2) 


*L—f>HUT-^Y'9JJ'bc, 

Here  we  adopted  the  regular  assumption  that: 

U  ~U. 

g*-) 

Therefore,  the  circulation  equation  (dimensionless)  is: 

V-~(Ub(U.i,-Ur) 

For  hinge  type  rotors,  we  assumed  that  the  blades  are  rigid,  its 
flapping  adjustment  coefficient  is  K  and  the  outward  displacement  is  e. 
For  hingeless  rotors,  it  can  be  assumed  that  the  blades  are  elastic  and 

it  is  treated  as  equivalent  quantities.  Thus 

+4, sin  i|> -Kfi. 

U,— 7+ l»sini> 

£/,*«—  A.+7+  Mcosiji -3,+  (r  —  e)*j^*  (3-3) 

and 

P«a~-y  i'-g~  2  (a.cosmJi+6.siniiTfr)^ 


(3-4) 


If  the  following  quantities  are  introduced: 


-K-jtr 


(3-5) 


a^,  bg- . . ,  we  get 


and  neglect  the  smaller  terms  containing  and  higher  than  second 

order  flapping  coefficients,  such  as 

f,.= — -Ma,  +  6?F +-|-Ma.  +  ) 

\ 

-^-Mo?  +  KaJ+ibs] 
ft#*-  ’  M6* —2atr  +  Kbtr^ 


(3-6) 


■,-J 
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It  is  worthwhile  noting  that  und^r  the  simplified  method  des¬ 
cribed  above  (that  is,  only  considering  the  induced  velocities  caused 
by  lower  order  and  same  order  terms  of  circulations),  the  derived  rela 
tion  between  induced  velocities  and  circulation  has  a  closed  form. 

By  solving  (2-5)  and  (3-6)  simultaneously,  the  various  orders  of 
harmonics  of  the  induced  velocity  can  be  obtained: 


r^E7.(«+4*»+-J-(»+4-i‘<«)) 

»<*+ l+rTWr5" 

\+AZ  l  3  3o  2  J  +  l  +  W?JU  +  .-iS  l 


A", /cos  a 


2a.r+Kblr  I 


+A«  (-f + 

H)  (3 


-A-,') 

10  y 


(3-7) 


where  4J,  ...are  defined  as  follows: 


.  2siti  a,  Aii-Al- 

’  1  +suia, 

*1 

„  2  +  2sin1a, 
"(l+sina,)* 

A\\-A\ 

„  2 cos  a, 

“i+sina, 

A\t  ~  A\, 

.  2— 2sin  a, 

0  1+sin  a, 

Aat^Ati 

0  2cos  a,  sina, 

'  (1+sin  a,)* 

AY, -A', 

2cos  a, 

#  (1+sin  a,)1  - 

AV.~AY. 

We  must  explain  here  that  in  the  deviation  of  induced  velocity  by 

integration,  if  infinity  appears,  then  it  is  necessary  to  replace  the 

lower  limit  by  r  instead  of  0.  F  is  the  relative  radius  of  the  blade 

o  o 

root  at  which  the  wing  cross  section  begins  to  show. 


Since  the  circulation  r  can  be  expanded  into  a  Fourier  series, 
the  airloads  can  also  be  written  in  the  form  of  a  Fourier  series: 

^--fw-^V+s  ((■§■).. «- (3-8) 

fit  *  l 

From  equations  (2-2),  (3-1)  and  (3-3) >  we  get 
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(M  --Kv-h*) 

““Hf”7  +_r1,f‘') 

(§“)„  r^-) 

(§i)„  — Hi‘“7  +4_‘‘f,‘ — r4*-) 

(3-9) 

Therefore,  the  relations  between  the  various  harmonics  of  the 
airload  of  the  blade  and  the  various  harmonics  of  the  circulation  are 
established. 

As  for  the  strain  coefficient  of  the  entire  rotor,  it  is  easy  to 
find  out: 

where  rQ  is  the  relative  radius  of  the  blade  root  at  which  the  wing 
cross  section  begins  to  take  shape;  r-^  is  the  loss  factor  of  the  blade 
tip  if  it  is  necessary  to  be  considered. 

IV.  FLAPPING  CONDITION 


Because  flapping  coefficients  are  included  in  the  expressions  of 
circulations  and  induced  velocities,  it  is  necessary  to  study  the  flapp¬ 
ing  motion. 

For  hinge  type  rotors,  according  to  [6],  the  flapping  motion  of 
the  blades  is  described  by  the  following  equation: 

.+0.Q*(/.+«5.)  -  (Af4).-  {M0),  (4-D 

where 

/.-  J*  (  r  -  0)tdmx  moment  of  inertia  of  a  single  blade 
e  around  the  horizontal  hinge 

S.-  f  *  ( r  —  * )dmt  moment  of  mass  of  a  single  blade  around 
J  *  the  horizontal  hinge 
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W«).-  J  r#  (  r  -  e  )dTl  strain  torque  of  a  single  blade  around 

the  horizontal  hinge 


gravity  torque  of  a  single  blade  around 
the  horizontal  hinge 
Written  in  a  dimensionless  form: 


It  can  be  further  written  as 


d*B  -  ~  _ 


(4-2) 


where 


l  + 


es 

(l-e)j 


S.-S./m,# 

(M  A)./m^Ri 
(Me) .  «■  (JWC)  ,/niiQ1!?1 


S“S./(l-e) 
Ma*~<Ma)./(  l-«) 
Me*  SS 


If  we  expand  MA 


then  we  get : 


into  a  Fourier  series 

aL-(m.),+  J  C (m^) cos »n|>  +  (Mj) « 

»-i 


6.(ir— v1)?*  (mJ« 


(4-3; 


(4-4) 


Therefore,  in  order  to  calculate  the  flapping  coefficients  aQ...an, 

b  .  we  must  first  solve  (M. )  - ,  (H. )  ,  (M. )  .  These  are  tedious 

n  a  o  a  no  A  nso 

integrations.  In  a  simplified  situation,  let-  us  assume  5  is  a  constant, 
e  »  0,  v1— l  and  note 

Y  —  ■  /j 

mt  2  / 

we  get 

“•-OTK-f+'T,l*)+A*  H" +-rl“)+-r« 

-  J » VW  -  J  j  ]  “T"  " 

fl‘  "  ["T"**1  +“J- A«M-  +  “  J  j  ~  J  o 

«-(J !  5-fW+K  -T  W'+~r°.»}/ (4-+-r>“) 

-  C«i JU±SrPt)/(~9l+  9») 

&*  -  (?.  A ~  <7t ~pj7 (?i + 9») 


where  j>t— - r  —  ^  5,.?W  +  -|-WS 

A-  -  J  -  J  *  - 7" 

•>‘-1 ^-*(-n+-rl“) 

— f1** 

Here,  we  noticed  that  amplitudes  of  the  higher  order  harmonics  of  the 
flapping  coefficients  are  smaller  than  those  of  lower  order  harmonics. 
In  determining  an  and  bn,  only  «►„  b„t,  <j^  ^ 
are  taken  into  account  and  a*„  b^u  a^t,  b^t 


In  addition,  we  get 


In  the  equations  of  the  flapping  coefficients  a„  a\,  6?—  and  strain 
coefficient  CT,  when  compared  with  the  classical  equations,  the 
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difference  is  that  terms  regarding  the  induced  velocity  integrals  are 
added.  Using  equation  (3-7),  we  carried  out  the  integration  and  fur¬ 
ther  obtained: 


r  Y  k  *  ,-lW  I  ,  1  ... _ d?'/m _ 4  A 

Ct  x~a-bld\3  (l+^S)  2  (l+.iil)  2  ( 1  +--1SK  i  +A\'t)  J 


(uli 


•^1,/cosa, 


l  7  (1+^«)  T  4  4  (  1  +  ^3)(  1  +A\*t) 

/ 

+x: 


-\LX 


) 


(  1 

A!,/cosa,  ..  ,Y| 

\  2  (1+^3) 

(1  +^8)(1  +Ali)  /J 

(4-8) 


„  vvf^  1  ,  1  _ ^/cosaj _  A 

o,  XY[fl\4  (i-MJ)  +4(  l  4  (  1 +.43)(  1 +^}{)  •  / 

■  ,  1  ».  -^tz/cosg) _  A 

+  \5(  1  +^3)  +  6(  1  +A\,,y  6  (  1  +v!3)(  1  +Ail)  J 

■  . J  1  ^t/cosa, 

+  \3(  1  +^3)  2  (  1  +  ^3)(  1  JrA\,t)  )\  3 

.  f  1  l _ 1 _ ^3,/cosa,  \ 

j  +Ay  +3(  !  +  ^}j)  3  (  l  +^3)(  i  +AH)  J 

(  1  ,  1  ^?,/cosa, _ ^ 

\4(  1  +A%)  +  4(  1  +A{',)  4  (  1  +^S)  (  1  +A[ J 


(4-9) 


_ f  1  ^lz/cosa, _ Y| 

Ul+^sr  2  (  1  +^J)(  1  -K4tf) 


^43,/coaa. 


r _ . _ 

UTT^iJ)  +  4  U  -r.43) (  1  -Mi*,)  4(  1  +.-13) 

6:-[3rrW)"a>+rrarnTT3ir)  (•••»«+•••«** 
+  o.2(M+-‘-^))]/(i(Ti3n5  +  4-] 


1  -u*  +  — *-4-]  I'  4—10) 

3  J 


(4-11) 


where  and  5. 


were  omitted. 


It  can  be  observed  from  equations  (4-8),  (4-9)  and  (4-11)  that 
under  the  same  flight  conditions  (*•»  **)  the  calculated  values  of 

CT  and  aQ  are  smaller  when  the  induced  velocity  variation  distribution 
is  considered  as  compared  to  those  obtained  with  a  constant  distribu¬ 
tion  of  induced  velocity.  However,  by  the  same  token  b3»  will  be  much 
larger.  Just  as  discussed  in  [7]  and  [9],  the  longitudinal  distribu¬ 
tion  of  the  induced  velocity  has  a  significant  effect  on  the  size  term 
b*  of  the  flapping  coefficient.  This  is  very  important  for  lateral 
control.  This  paper  for  the  first  time  introduced  the  analytical 


expressions  of  CT  and  flapping  coefficients  which  are  only  dependent 
on  the  flight  parameters  but  has  taken  induced  velocity  variation  dis¬ 
tribution  into  account. 


In  addition,  we  have: 

a,-  (QttPn  +9i  J>**)/<9«>9«  +  9i*9t,)  (  4_12  ) 

where  the  definitions  of  Put  Pw  9u»  9u»  9n»  9»*  are  shown  in  [10]. 


V .  AIRLOADS 

Substituting  the  expression  of  harmonics  of  the  induced  velocity 
(3-7)  into  that  of  the  circulations  (3-6),  we  can  get  the  following 
matrix  form  circulation  equation:  r  ji*  'i 


(5-1) 


where  the  elements  of  the  matrix  [Q]  are  shown  in  Appendix  1  in  [10]. 


Then  substituting  circulation  equation  (5-1)  into  the  airload  equa 
tion  (3-9)  for  blades,  we  finally  get 


where  the  matrix  [P]  is: 


VI .  EXAMPLE 


In  order  to  explain  this  method,  we  used  the  Yen-2  helicopter  rotor 
blade  as  an  example  to  calculate  the  flapping  coefficients  under  the 
seven  flight  conditions  of  i»-o.05,  0.075,  0.10,  0.125,  0-15,- 20,  0.24 
and  the  airload  where  y  ®  0.20. 


R- 5  «n 

AO- -0.1396  rad  ' 
if —  0.3 

m,  =  2. 755  kg-secVm 
P«— 0.108  kg- sec*/ m 


5—0.0486 
e  — 0.014 
h  -  3 

Q—  37.48  rad/sec 
o.-5. 73 


The  flight  parameters  were  obtained  from  the  balance  calculation 
of  the  craft.  Using  y  *  0.20  as  an  example 

V,— 0.2053  cosa,-  0.9741 

OJ— 0.2409  rad  XJ- -0.02494 


Based  on  equations  (4-9),  (4-10),  (4-11),  (4-12)  and  (4-13)  to 
calculate  the  flapping  coefficients,  we  obtained  the  curves  showing 
the  relations  between  aQ,  aj,  bj,  a2,  b2  and  the  advance  ratio  y  as 
Figures  1,  2,  3,  4  and  5  respectively. 


For  comparison  purposes,  the  flapping  coefficients  obtained  with 
constant  induced  velocities  are  also  shown  in  these  figures.  We  can 
see  that,  as  we  described  before,  the  aQ  curve  with  induced  velocity 
variation  distribution  is  lower  than  that  with  constant  induced  velocity, 
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The  aj  curves  for  both  types  of  distribution  are  almost  identical.  How¬ 
ever,  the  b|  curves  of  both  distributions  are  quite  different.  The 
former  is  greater  than  the  latter  especially  exhibiting  peak  values 
in  the  low  velocity  region.  This  phenomenon  has  been  observed  many 
times  in  experiments.  a2  and  (-b2)  curves  are  similar  to  those  of  b£ 
but  they  are  less  than  l/10th  of  that  of  bj.  Prom  these  observations, 
it  is  reasonable  to  neglect  higher  order  flapping  coefficients  in  the 
calculation  of  lower  order  flapping  coefficients. 


Finally,  based  on  the  equations  of  induced  velocities  and  various 
orders  of  harmonics  of  blade  airloads  (3-7)  and  (5-2),  the  distributions 
Of  5„  5...  8,„  5„,  and  ,  (*£■)_.  ($■)„'  ’  (€*-)„ 


>u 

The 


along  the  radius  under  the  condition  u  *  0.2  were  calculated.  The 
results  are  as  shown  in  Figures  6-10  and  11-15.  In  order  to  verify 
the  validity  of  this  simplified  method  in  this  paper,  the  varous  har¬ 
monics  of  airloads  of  rotor  blade  obtained  using  numerical  integration 
[11]  are  also  shown  in  Figures  11-15.  It  is  clear  that  the  results  of 


airloads  obtained  using  two  different  methods  are  the  same.  In  addition. 
Figure  16  shows  the  curve  of  variation  of  airload  along  the  azimuth 
angle  at  different  radial  distances.  The  trend  of  these  curves  is  also 


very  similar  to  that  discussed  in  [1] 


Fif.l  Coaist  ancle  °s  advance  ratio  “ 

1 —  equation  in  this  paper 

2 —  classical  equations 
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Fig. 2  Flapping  coefficient  °*  ts.  advance  ratio  u 

1 —  equations  in  this  paper 

2 —  classical  equations 


Fit .  3  Flapping  coefficient  v*.  advance  ratio  t* 

1 — equations  in  this  paper;  2 — classical  equations 


Fig.  4  Flapping  coefficient  «» 


advance  ratio  F 

1 —  equations  in  this  paper 

2 —  classical  equations 


Fig.  S  Flapping  coefficient  *»  **•  advance  ratio  I4 

1 — equations  in  this  paper;  2 — classical  equations 


Fig.  8  Distribution  of  th«  induced  Telocity  along 
radius  for  I*  >0.20 


Fig. 7  Distribution  of  the  induced  Telocity 
»u  along  radius  for  I1  >0.20 


Fig. 4  Distribution  of  the  induced  Telocity  *\t  along  radius  for  »  *  0  20 


Fjg.  9  Distribution  of  the  induced  velocity 
*u  along  radio*  for  I*  "0.20 


Fig.  10  Distribution  of  tbe  induced  velocity 
*u  along  radius  for  I1  *0.20 


Fig. 11  Distribution  of  tbe  blade  airload  ( along  radius  for  I1  *0.20 

1 — equations  in  this  paper;  2 — numerical  solutions 


Fig.  12  Distribution  of  tie  blade  airload^  ^  )  along  radios  for  l*  *0.20 

1 — equations  in  this  paper;  2 — numerical  solutions 


Fig.  14  Distribution  o(  tbe  blade  airload^  -j~-  along  radius  for  (**02 

1— equations  in  this  paper;  2 — nunerical  solutions 


VII  CONCLUSIONS 


The  major  conclusions  of  this  work  are  summarized  in  the  following: 


(1)  Based  on  the  classical  vortex  theory  of  rotor  and  the  blade 
element  theory,  a  closed  form  relation  between  induced  velocity  and 
circulation  has  been  established. 


(2)  For  the  first  time  analytical  expressions  for  the  flapping 
coefficient  and  blade  airload  which  not  only  take  effect  of  the  variation 
distribution  of  the  induced  velocity  into  account  but  also  only  depend 
on  blade  parameters  and  flight  parameters. 


(3)  The  method  to  estimate  rotor  airload  introduced  in  this  paper 

is  simple  and  suitable  for  engineering  applications. 
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A  SIMPLIFIED  METHOD  FOR  PREDICTING 
ROTOR  BLADE  AIRLOADS 

Wang  Shieun  and  Xu  Zhi 
( Nanjing  Aeronautical  Institute) 

Abstract 

At  present,  a  simplified  approach  to  the  prediction  of  rotor  blade  airloads 
is  urged  to  be  developed  in  the  engineering  application. 

In  this  paper,  firstly,  relations  of  first  two  harmonic  induced  velocities  to 
the  lower  and  same-order  harmonic  circulations  are  obtained  from  the  generalized 
classical  vortex  theory  of  the  rotor.  Then,  based  on  the  blade  element  theory, 
a  closed  form  of  equations  for  circulation  is  established  and,  by  taking  the 
flapping  condition  into  account,  simplified  formulae  for  predicting  rotor  blade 
airloads  are  set  up.  In  particular,  expressions  of  flapping  coefficients  are 
derived,  including  the  effect  of  variable  induced  velocity  distribution  but  in 
terms  of  blade  parameters  and  flight  parameters  only. 

Finally,  a  calculation  of  a  typical  example  is  made  and  comparisions  of 
airloads  with  those  from  the  more  accurate  numerical  solution  are  shown  that 
the  present  method  is  fairly  suitable  for  aerodynamic  analysis  and  preliminary 
design  of  helicopters. 
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STUDY  ON  PRESSURE  DISTRIBUTION  ON  ROTOR 

BLADES  WITH  THREE-DIMENSIONAL  NON-STEADY  THEORY  OP 
COMPRESSED  FLUID 

Jiangxi  Aeronautical  Society,  LI  Zenhao  and  Ruan  Tlcusen 
ABSTRACT 

A  solution  to  the  pressure  distribution  on  rotor  blades 
tnree-dimensional  nonsteady  theory  of  compressible  fluid  is 
presented  in  this  paper  in  the  case  of  continuous  wake-surface 
and  instant  speed  forward  motion  of  a  helicopter.  In  order 
to  satisfy  the  wake  conditions,  it  is  assumed  that  the  press¬ 
ure  doublets  move  along  the  wake  surface  instead  of  along  the 
actual  tracks  of  blades.  By  adding  the  moving  pressure  doub¬ 
lets,  an  integral  equation  of  the  three-dimensional  nonsteady 
compressible  fluid  with  superior  singularity  is  obtained  when 
the  blades  are  in  complex  motion.  A  method  to  treat  this 
superior  singularity  was  obtained  to  make  it  a  continuous  inte¬ 
gral  equation  with  super  extreme  values.  Furthermore,  the 
numerical  solution  of  this  equation  was  also  obtained.  The 
accuracy  of  the  method  and  the  equation  has  been  proven  based 
on  simple  examples  using  a  small  computer. 

I.  THE  PROBLEM 

Up  until  now,  the  calculation  of  helicopter  rotor  loads  by  tradi¬ 
tion  adopts  a  vortex  line  simulated  wake  of  Irregular  or  non-singular 
shape  and  uses  the  nonsteady  incompressible  fluid  Pierre-Savor  equation 
to  obtain  the  induced  velocity  at  the  blades  and  determines  the  lateral 
distribution  based  on  classical  blade  element  theory.  The  calculated 
results  usually  show  a  lower  total  lift  at  higher  speeds  or  a  non-zero 
load  at  the  root  ahd  tip  of  the  blade.  Even  in  some  cases,  it  is 
necessary  to  manually  shift  the  wake  position  to  approach  the  experimen¬ 
tal  results.  The  error  for  higher  order  load  is  also  large.  Higher 
order  load  is  the  source  of  helicopter  vibration.  Therefore,  it  is 
Important.  The  reasons  for  the  above  discrepancies  lay,  of  course,  in 
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the  assumptions  made.  It  is  not  uncommon  to  see  the  use  of  acceleration 
potential  to  solve  nonsteady  aerodynamic  problems  for  fixed  wing  air¬ 
craft  but  it  is  rarely  used  in  treating  rotors.  With  regard  to  [1], 
we  will  make  some  comments  in  the  following.  As  for  [2],  it  only 
carried  out  some  exploratory  discussion.  This  paper  intends  to  over¬ 
come  the  inadequacy  in  the  area  and  to  develop  a  new  rigorous  accelera¬ 
tion  potential  method  with  nonsteady  theory  of  compressible  fluid  taken 
into  account  to  calculate  the  airloads  of  rotors. 

Fixed  coordinate  systems  on  the  ground  o-x,  y,  z  and  o-r,  0,  z  are 
chosen  with  z-axis  positive  in  the  upward  direction.  At  time  t  =  0, 
the  center  of  the  blade  shell  coincides  with  the  original.  The  blade 
on  which  the  load  is  calculated  is  located  downwind  (<J>  *  0).  x-axis  is 
in  the  opposite  direction  to  U^.  is  the  constant  forward  speed  of 

the  helicopter.  Assuming  that  no  shock  wave  and  no  separation  exists  in 
the  flow  field,  the  blade  only  undergoes  a  flapping  elastic  motion  as 
shown  in  Figure  1. 


Because  the  attack  angle  curvature  and  thickness  are  not  large  when 
the  blades  are  working,  the  disturbance  generated  due  to  the  motion  of 
the  blades  is  very  minute.  At  this  time  the  acceleration  position 
<j> . — exists,  and  the  following  equation  is  satisfied 


**<f>  ,  ^4>  .  **4> _ l_  n 

dx*  ay1  +  iz*  a *  a;1 


(1) 


The  following  is  to  determine  the  relation  between  the  acceleration 
potential  and  velocity  potential  $.  According  to  nonsteady  Bernoulli 
equation,  the  acceleration  potential  at  reposity  is  zero. 


received  in  May  1981 


but  because: 


dj  _  *5 

dt  it 


+o. 


^  i  ^ 

— -  —  -r  g.  — ; — +»,— - — 

iX  *  9y  iz 


i£ 

at 
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therefore 


When  a  helicopter  is  suspending,  the  axial  induced  velocity  is 
the  largest.  Based  on  the  momentum  principle,  (T/2npR*)l'*t  where  T 

is  the  weight  of  the  flying  object.  Therefore,  7*/(4P«i?*)*' 

but  S/(4ttR^)  in  reality  is  between  0.07-0.1. 


T  S  bp  S_ 

ipst:  u*&r~  p  u*#) 


Therefore,  o*/2  is  far  less  than  (P-pm)/Pm.  So,  under  small  disturbance 
conditions : 

A  — - j^-3— t. 
v  *t  dt 

or 

**  J-l  (2) 

only  t  is  a  variable  in  the  integration,  x,  y  and  z  do  not  vary. 


It  is  clear  from  the  above  discussion  that 
are  not  identical  when  »V2  cannot  be  neglected. 


di 


dt 


and  grad  4  ~~ 

at 


When  an  acceleration  potential  method  is  used,  it  must  satisfy 

(1)  at  infinity,  the  pressure  is  p^,  and  <p  =  0 

(2)  the  acceleration  potential  4>  satisfies  the  wave  equation  (1) 

(3)  the  velocity  potential  satisfies  equation  (2) 

(4)  the  wake  created  by  the  motion  of  the  rotor  produces  a  wake- 
surface  in  the  flow  field.  The  wake-surface  is  dragging  out 

of  the  rear  fringe  of  the  blade.  At  any  instance  across  the 
wake-surface,  the  velocity  potential  i  is  non-continuous  with 
a  sudden  jump  in  quantity 

(5)  on  the  surface  of  a  moving  blade,  boundary  conditions  are 
satisfied.  This  means  that  the  normal  component  of  the  flow 
velocity  at  any  point  on  the  surface  at  any  instance  is  equal 
to  the  instantaneous  velocity  component  at  this  point. 

II.  MOVING  DOUBLET  ACCELERATION  POTENTIAL  INTEGRAL  EQUATION 


Doublets  are  distributed  in  the  flow  filled  to  satisfy  the  five 
conditions  mentioned  before.  Doublets  are  expressed  as  (C,n,S)  and 
the  calculation  points  in  the  flow  filled  are  (x,y,z).  In  order  to 
satisfy  condition  (4),  doubtlets  move  along  the  wake-surface.  Let  us 
assume  at  timet^,  the  doublet  is  located  at  R1  and  at  t2  it  moves  to 
R2,  then  4>  has  a  discontinuity  across  the  surface  from  R1  at  t^  to  R2 
at  t2*  However,  <t>  does  rot  jump  at  R^  when  time  is  t2.  As  for  f, 


24 


from  equation  (2),  it  is  clear  that  $  is  continuous  because  <J>  is  con¬ 
tinuous  at  any  point  outside  the  wake-surface.  When  crossing  the  wake 
surface,  $  also  leaps.  This  leap  occurs  at  between  t^  and  tg*  There¬ 
fore,  moving  the  doublet  along  the  wake  surface  will  make  the  motion 
satisfy  condition  (4). 


The  moving  doublet  acceleration  potential  is 

/(Qa, 
<i.a  Qm 


1 


*(*,  y,  t  )“ — ~ 


a 


(3) 


nQ  is  the  axial  direction  of  the  doublet,  from  junction  pointing  toward 
the  source.  In  this  paper,  the  junction  is  placed  on  the  wake  surface 
or  under  the  blade.  Partial  derivatives  are  made  with  respect  to  (£,n» 
C).  f(x) — doublet  strength  function;  a^ — undisturbed  sonic  velocity. 

a 


t  -■ 


a — distance  between  doublet  point  and  flow  field  point 


a«>  c *  —  £(t)]*+Cy-n(T):*+C*“t(^)3I 

sn'c-o-cy-nmo-c*- 


(4) 

(5) 

(6) 


After  substituting  (3)  into  the  wave  equation  and  repeatedly  carrying 
out  mathematical  operations,  it  is  possible  to  prove  that  (3)  satisfies 
the  wave  equation  C,n,S  are  the  coordinates  of  the  doublet  when  it  is 
moving  along  an  arbitrary  curve.  Prom  (3)— (6),  it  is  demonstrated  that 
the  parameter  a^  which  represents  the  effect  of  compressibility  is 
incorporated  into  the  pressure  verturbation  in  a  complicated  form.  It 
does  not  follow  a  correction  such  as  v'\  —  Mi,  as  used  in  Planto's 
correction. 


Let  the  number  of  blades  be  MM.  The  potential  difference 
between  the  mth  blade  and  the  blade  on  which  the  load  is  calculated 
is  tpm,  i£m  *  (m-1)  2ir/MM.  Then  introduce  the  blade  coordinates,  p,  0Q; 
r,  ij^0  as  shown  in  Figure  1.  0Q,  (i^q)  is  positive  when  pointing  toward 
the  front  fringe.  The  doublet  is  moving  on  the  surface  of  the  blade 
when  traveling  from  the  corresponding  position  on  the  blade  surface  to 
the  rear  fringe.  For  simplification,  let  us  assume  that  it  moves  on 
the  z  ■  0  plane.  Beyond  the  rear  fringe,  the  motion  of  the  doublet  is 
affected  by  the  perpendicularly  downward  wash  velocity  w  and  it  is 
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Fig.l  Blade  tad  coordinate 


Fi(.2  Blade*  positions  at  different  moment 


dragged  downward.  Assuming  w  is  a  constant,  it  corresponds  to  a  non 
Irregular  wake. 


Let  the  time  for  the  doublet  to  move  from  the  doublet  point  to 
the  rear  fringe  be  tQ1, 

f,a-(  Psin  8.+0.75  b  )/CPQ+l/-sm(<*>  +$.)}  ( 7 ) 

According  to  Figures  1  and  2,  we  get 

*  "  r  C08(cp  +’!>»)*  y  *  r  sin(<p  +$,),  (8) 


4-- 


C/-t+Pcos(QT+^+e,),  Ti-Psin(QT+^+0,),  t  -  »[r  /,,)]  (9) 

0 


<p  w 

~cT>  x  o 


constant 


(10) 


o  *•» 


Ui>  *  >-oo 


°  -  jf  rcos(q>  +^,)-  Pcos(Qt+*..+8,) +£/.(% 

+  C  r  sin(«P  +  *,)-  Psuj(Qt+^..  +  04))*  +  (  z  -  w  C  t 

0--  -  [  r  cos(  <P  +  *,)  -  P  cos  (Qt+^+9.)  +uj(  T - ^-)jc-t/.-PQsiii(QT+i. 

+ 0.)  ]  -  C  r  sin (  <P  + *,)  -  P  sin  (Qt + + 0t) ) PQcos(Qt + $_+ 9#) 


(ID 

(12) 


~(~5 — '“)]}■ 

a 

In  order  to  find  the  partial  derivatives  in  the  normal 

direction,  it  is  necessary  to  obtain  the  normal  directional  number  of 
the  slanted  helieoidal  surface  forward  by  the  trajectory  of  the  doublet 
Solving  (9)  simultaneously,  we  get 


I  +  — *..)+  n  +  — *#i)+v-+9.]  d3) 

Talcing  partial  derivatives  of  (13)  with  respect  to  (,  n  and  c,  we  get 
the  normal  directional  numbers 


u>sin(Qr+ —  a»coa(Qr  +  ii.  +  9,),  PQ  +  C/.sin(QT  +  1}„,4-0t) 


(14) 


After  that,  we  must  use  lim  --jj-  /•  ,  and  carry  out  compli¬ 

cated  operations  to  get  the  following: 

..  *  /  (  *  )<*»  ,  (  ,  x  l*  +  CQA-bc)/(.bd-  g) 

,!T5  dz  **•  *  '  '  **~bd-y/vr+bt* 

r  ( t )  QA  bf+b«-  *Sbd  ./'(  T)  (15) 

a.  0*53  v/u>1  +  b*1  a*  aVta>/ur  +  6«?1 

where  b«=-PQ+t7-s”>(QT+i-+0») 

—  (  jc  —  | )  wsuUQt-*- 9,)  +■(>  —  n)(— u»  )cos(Qt-  v„  +  9,)  —  i 


a,  are  equations  (11)  and  (12)  where  z  «  0. 


When  the  doublet  is  moving  arbitrarily,  the  relation  between  the 
pressure  differential  and  the  acceleration  potential  strength  is 

-£z^--/(t)  (17) 

when  the  helicopter  is  flying  forward  at  a  constant  speed  U^,  P^own“ 
Pup  «  Ap  can  be  written  as:  (Ii*  *ii  are  the  dimensionless  orthogonal 
coordinates  of  the  blades). 

A  P  -  2  ’ii)*'*'0’4**’ 


(18) 


Therefore,  the  velocity  potential  generated  by  the  doublet  is 

1  f*/Q  a  ,, 

♦“ISptJ-oo  2  a  +QJo-d1 

»«0 

The  above  equation  is  for  a  corresponding  infinitestimal  area.  For  the 
entire  rotor,  doublets  are  distributed  on  all  the  blade  surfaces. 
Therefore,  we  have  to  integrate  with  respect  to  the  surface  of  the 
blade  and  to  add  the  number  of  blades  to  it.  We  get 


4  4*p.  J-J-1^2  A*,(**» 

*  ■  0 

^  (9/Q  a  expC»»(QT  +  ^.)]  Jj 

mm  l  J  —  “  +Q-/I- 


(19) 


where  S  is  the  area  of  a  single  blade, 


By  taking  the  partial  derivative  on  both  sides  of  equation  (19) 
with  respect  to  z  and  using  condition  (5)  together  with  assuming  that 
the  normal  velocity  of  the  arc  surface  on  the  blade  is  V  ,  we  get 

*2  vfj'Z  -sr 

When  the  integration  is  carried  out,  the  relation  between  rj,  and  P,  0, 
must  be  used.  Based  on  Figure  3,  we  get 


— 2  P  sin  0B 


a  r  =  m  2  P  COS0,  .  ZRS 

0.5,  'll  *  j  1  1  £ 


P,*0.5  y/ Lz^ti,  +  1  +  +  6,(|1  +  0.5)»,  0o  =  tg"  iJ~  +  i  +  2RS\  | 


(21) 


III.  THE  TREATMENT  OP  ADDITIVE  SINGULARITIES  OF  WAKES 


Equation  (20)  is  a  singularity  integral  equation  of  rotor  load 
under  nonspeady  compressible  conditions.  Its  core  is  an  integral 
which  has  not  yet  been  solved  to  date.  The  double  integral  on  the 
outside  and  the  partial  derivatives  for  finding  out  the  extremes  are 
also  very  difficult.  In  addition,  pressure  differential  between  the 
front  and  rear  fringe  and  boundary  conditions  at  the  tips  and  roots 
of  the  blades  must  be  satisfied.  In  order  to  find  a  solution  for  (20), 
we  have  to  study  the  experience  of  earlier  work  in  the  research  of  air¬ 
craft  wings. 

In  [1],  the  method  used  to  solve  the  similarity  equation  s  is  to 
use  z  -*■  zQ  to  substitute  z  -*■  0.  zq  is  not  equal  to  zero  and  it  is  a 
small  quantity.  Based  on  the  calculated  results  in  that  paper,  it  is 
0.035  times  the  arc  length.  Thus,  a  continuous  integral  equation  is 
obtained.  That  paper  still  used  a  numerical  method  to  find  the  deriva¬ 
tives  by  choosing  a  step  length  of  0.005  times  the  arc  length.  It 
appeared  to  be  quite  reasonable  significally  because  the  difference  in 
normal  direction  velocity  V  at  z  ■  0  and  z  »  0.035  arc  length  is  not 
too  large.  However,  this  arrangement  in  practice  has  a  significant 
amount  of  artificial  arbitrariness.  Using  the  simplest  example  of  the 
incompressible  zero  order  wing  load  problem  (in  the  analysis  of  this 
method,  there  is  no  practical  difference  between  wings  and  rotors), 
because  of  the  singularity  of  the  integral,  it  is  not  necessary  to 
numerically  find  the  derivatives.  It  is  possible  to  move  the  calcula¬ 
tion  of  limits  and  partial  derivatives  into  the  integral.  Thus,  the 
kernel  becomes  After  the  calculation,  it  becomes  l-t**/**  where 
r  is  the  distance.  When  the  coordinates  x  and  y  of  the  flow  field 
point  and  the  doublet  are  identical,  which  means  that  we  are  looking 

for  the  point-to-point  effect,  r  -  z  .  For  the  kernel  of  the  integral, 
3  0 

it  becomes  -2/Zqq.  Therefore,  the  main  diagonal  elements  of  the  matrix 
completely  depends  on  zQ.  When  zQ  changes  from  a  small  number  to  a 
large  number,  the  airload  can  change  from  very  large  to  very  small. 
Thus,  for  different  rotor  and  flight  conditions,  it  is  necessary  to 
obtain  the  load  data  before  zQ  can  be  decided.  The  meaning  of  solving 
the  unsteady  compressible  rotor  load  problems  is  completely  lost. 


In  the  nonsteady  compressible  wing  problems,  the  lim-v-  can  be 

j»0  ” 

moved  into  the  double  integral  and  the  kernal  function  becomes 

.  f  **  «xpCio(X  -A/v'x*+B»y+3V  )) 

i-o  *■  J  —  77+^v+SV - x 


After  integrating  this  equation  and  isolating  the  singular  point, 
then  it  is  possible  to  find  the  Hadamard  principal  value  of  the 
double  integral.  If  we  can  separate  (20)  into  two  parts,  one  with 
singularity  and  one  without  it  and  make  the  part  with  singularity  to 
have  the  same  form  as  the  one  for  wings,  then  we  can  use  the  results 
of  wings  to  make  the  rotor  loads  as  the  wing  loads  plus  some  corrections 
This  is  the  method  used  in  this  paper. 


The  following  is  a  derivation  of  the  solution  to  this  equation. 
Considering  that,  among  the  contributions  of  all  the  blades,  the  effect 
of  other  blades  on  the  one  being  calculated  is  much  smaller  due  to  the 
distance.  It  can  be  assumed  tabt  the  wakes  of  other  blades  fall  off 
from  the  doublets  which  correspond  to  m  /  1  at  tQl  =*  0.  Equation  (20) 
can  be  written  ae 

V.(r,  *t ,  cP)-~-  4Jtp^  lim  ~  J  _  j  J  _  j  J  A/>„(I„  n,) 

n  »  0 

xjjf  \*/a  j-  -i.  eim"  dt  (22) 

f<*/Q  d  e"*'  ) 

J  T/Q-f„,  (nt«l)  a  J  *  * 

V\x\  V™  are  continuous  and  without  singularity.  Let  us  move  lim — 

(3)  ‘ 

into  the  inner  integral.  contained  all  the  singularities  and 


transform  ----- 


—  f  «-f'. 


4  4JTP.  iLm0  *z  J-i  S  ^l) 

*  v  n  »  0 

f  0  wptwCP+OO) 

J  -i„  K  a  +0 m/a- 


(23) 
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Let  us  transform  r,  i|>Q;  p,  6Q  into  the  blade  coordinates  x^,  y^, 
5^,  n^,  according  to  equation  (21) 


P  cos  9 


Pnne,--|» - j“.  rani.--*, - j- 

4 

W+J^+RS,.  r  cmibt-ylJr~-+RS 


(24) 


Furthermore,  because  in  the  -t  ^  to  0  integration  limits  the 
minimum  value  of  ftt  is  around  0.04  (the  blade  has  a  high  span  to  arc 
ratio;  usually  it  can  reach  16),  thus  it  is  possible  to  let 

0*T*  „ 

«  ••  A*  A  -r  T  4-  4  trt  n  1  o  a  ^  /s  atm  4  4“  T  7  ffl 


COsQt—  i  —  - 


sin  Qt=»Qt. 


It  is  also  possible  to  omit  Um cos<=p 


and  all  the  higher  terms  than  t2  to  obtain 

o  -  ✓  U,  - 6, + 1/, T)* • +  (v , — Tl, y '  +  (  r  -  t  )*" 
0L-i;1(*1-t1+i/,T) 

L^  —  pQ+LLsin  tP 


(25) 

(26) 
(27) 


Under  these  assumptions,  it  is  possible  to  express  t  explicitly 
with  t.  Its  method  is  similar  to  that  of  the  wing.  Finally,  we  get 

r  -557-  f-! 


where  3 


\  v  f°  f.iQf,  ,  («,~l,)C/, _ F  V  dt 

1  —AT*,  -1.,  eipi'  P‘l<+ 

^-^(*,-6,+y,0,  +  3I(y,-Ti,),  +  3,(  *  -  t  )• 


(28) 


(29) 


The  symbol  represents  the  need  to  obtain  the  Hadamard 

principal  value  of  the  double  integral  because  the  operations  of  cal¬ 
culating  the  derivatives  and  limits  are  moved  inside  the  double  integral 
sign.  Let  us  write  the  inner  integral  as  J  and  make  £,+(/,/-  X, 

„pr.  «KS-t,n 

1 - 177"“ - “31?'  (3°  • 

f„_tl  ]} 

K  i/x* + 3*(y  1 — Tii)1 + 3V 

The  form  of  integral  s  in  (30)  has  been  solved  in  [4],  The  inte¬ 
gration  limits  are  from  -»  to  xQ.  Using  equation  (D8)  in  that  paper 
and  making 


/(X)rfX- /(X),X-J-—  /(X)rfX  (31) 


this  corresponds  to  the  tangential  motion  along  the  blade  added  or 
subtracted  by  a  term  which  is  the  wake  from  to  the  rear  fringe  of 
the  blade.  In  addition,  by  noticing  that  the  terms  irrelevant  to  the 
upper  limit  of  the  integral  in  (D8)  are  cancelled  in  the  subtraction 
we  get 

•VC.  *>-- f- 
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Therefore, 
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where 
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IV.  METHOD  OP  CALCULATION 


Due  to  the  simplicity  of  the  blade  surface,  a  series  method  can 
be  used  which  has  a  lower  demand  on  the  memory  of  the  computer.  Let 

Ap«(5i»  *ii)! 


—■ J  ^:v,( i;) 


v  r+T 


p » i  ?  ■  i 


(35) 


where 


5-^t-  1.  gA x-v  i  4-IJ, 


V'  1 


(36) 


After  the  substitution  of  (33)  and  undergoing  a  dimensionless  process, 
a  set  of  linear  algebraic  eauations  regarding  u„.t  is  obtained. 

In  the  coefficient  matrix,  the  part  containing  the  singularity  is  ob¬ 
tained  using  the  method  described  in  [4].  The  finite  integral  in  that 
equation  can  be  approximated  by  (21)  in  [4].  But  the  calculation 
involves  complex  numbers.  In  this  case  the  complex  expression  of  a 
Fourier  series  to  replace  (18) 


A  P» 


J  ~~  (( Ap«,i  Ap«,,)e"<0,***>  +  ( Ap«.i  +»’Ap«,i)c”",0,**"J) 


(37) 


The  calculation  has  to  be  carried  out  only  with  respect  to  the 
term  corresponding  to  .  This  is  because  the  V  value  obtained 

from  the  calculation  of  the  other  term  must  be  a  complex  conjugate. 
Thus,  if  a  certain  order  is  discussed,  then  we  get 


F,-(o  +*6)— (Ap.,i  — «Ap.,t)  +  ( o  - ib )  -y-(A/>.„ +iA/0 
“  0  Ap.„  +6Ap.,i 


(38) 


where  a  and  b  represent  the  real  and  imaginary  parts  of  the  singularity 
portion  respectively.  Therefore,  only  by  using  the  term,  the 

real  part  calculated  corresponds  to  the  cosine  component  and  the  imag- 
nary  part  corresponds  to  the  sine  component. 


Equation  (3*0  is  a  non-singular  portion  which  is  a  continuous 
function.  However,  its  calculation  is  very  difficult  because: 

Firstly,  since  for  the  factors  containing  eima'  ft  reaches  about 
20,  therefore,  the  inner  integral  is  a  large  parameter  oscillating  type 
of  integral.  When  n  increase,  the  value  of  nft  increases  rapidly  which 
makes  the  number  of  integral  knot  point  to  increase  quickly.  The  number 
of  integrals  also  increases  rapidly  which  makes  it  impossible  to  cal¬ 
culate  higher  order  load  in  practice.  In  order  to  solve  this  problem, 
a  simple  function  is  used  as  an  extrapolated  function  for  the  function 
to  be  integrated  at  n  *  0  to  obtain  the  extrapolation  coefficient. 
Finally,  in  the  extrapolation  region  the  integral  is  preciously  obtained 
by  using  the  factor  which  is  the  product  of  three  sample  functions  times 
eimv  .  This  method  makes  it  relevant  between  the  calculation  of  the 
function  value  and  load  order  n.  It  is  then  possible  to  obtain  the 
highest  order  of  load  without  increasing  the  computer  time. 

Secondly,  the  function  to  be  integrated  has  very  strong  peak  values 
at  «-l.  y-Ti,  i  When  m  /  1  it  also  has  a  strong  peak  but  the 

position  is  hard  to  determine.  In  order  to  ensure  that  the  calculated 
value  is  not  covered  by  errors,  we  divided  into  subregions  between  [-1, 

1]  for  the  integral  in  the  span  and  used  y  =  n  as  one  end  of  an  interval. 
Then,  Gauss-LeChardin  model  is  used  to  get  the  solution. 

In  compiling  the  program,  it  is  necessary  to  minimize  the  work 
amount.  If  the  8th  order  load  on  the  lift  surface  of  a  four  blade 
system  is  to  be  determined,  the  number  of  integrals  can  reach  over 
100,000.  In  reality,  it  is  possible  to  eliminate  the  number  of  arc  and 
span  series  from  the  number  of  integrals  and  to  eliminate  the  load  order 
number  based  on  the  extrapolation  method  mentioned  above,  the  number  of 
single  intergrals  can  be  reduced  to  about  5000. 

This  method  can  also  be  solved  using  a  grid  method. 

The  azimuth  angle  <j>  is  to  provide  convenience  for  Fourier  analysis. 
After  performing  Fourier  analysis  on  both  the  matrix  and  V  ,  we  can  get 

a  ( 39 ) 
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Here  the  elastic  flapping  motion  of  the  blade  is  taken  into  considera¬ 
tion.  V  must  be  solved  from  the  differential  equation  in  the  dynamics 
z 

of  the  elastic  flapping  of  the  blade.  Thus,  (39)  and  this  equation 
are  coupled  as  frequency-state  for  a  solution. 

V.  SIMPLE  EXAMPLE 

The  original  data  of  this  example  is  taken  from  [5]  Table  15  which 
represents  the  flight  test  of  the  American  H-3**  helicopter.  Because 
of  the  inability  to  obtain  the  necessary  funding,  only  the  Planto- 
Glowatt  corrections  for  incompressible  working  conditions  are  computed. 
The  number  of  calculation  point  was  also  minimized.  It  is  carried  out 
on  the  slowest  model  121  computer  available  in  the  country.  Only  four 
azimuth  angles  were  taken  (therefore,  only  first  order  harmonic  analysis 
was  made  on  the  load).  As  for  the  lift  line  theory,  in  order  to 
reduce  the  work  load  for  the  computer,  we  increased  the  length  of  the 
portion  of  blade  without  any  wing  shape.  Thus,  the  load  at  the  root  of 
the  blade  has  a  higher  error.  The  results  are  shown  in  Figure  4. 

Based  on  the  comparison  between  theory  and  experiment,  the  zeroth 
order  total  lift  and  distribution  coincide  very  well  with  experimental 
results.  Under  the  conditions  that  the  number  of  azimuth  angles  is 
few  and  the  calculation  points  are  coarse,  the  first  order  load  ob¬ 
tained  agrees  very  well  with  the  experimental  results. 

The  major  problers  in  this  paper  were  guided  by  Professor  Low  Shi- 
chuin  of  Northwest  Industrial  University. 
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Figure  4.  The  comparison  of  theory  result  with  measured  in  flight 

1 — aerodynamic  lift  per  unit  span  (kg/m);  2 — theory;  3 — experimental 
4 — blade  span;  5 — (a)  zeroth  order  load;  6 — (b)  first  order  sine; 

7 — first  order  cosine 
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STUDY  ON  PRESSURE  DISTRIBUTION  ON  ROTOR  BLADES 
WITH  THREE-DIMENSIONAL  NONSTEADY  THEORY 
OF  COMPRESSIBLE  FLUID 

Li  Zhenhao  and  Ruan  Tianen 
( Jiangxi  Aeronautical  Society ) 

* 

Abstract 

A  calculation  of  pressure  distribution  on  rotor  blades  with  three-dimen¬ 
sional  nonsteady  theory  of  compressible  fluid  is  presented  in  the  case  of 
continuous  wake-surface  and  forward  motion  of  a  helicopter  at  a  constant 
speed.  An  acceleration  potential  equation  is  derived.  A  fundamental  solu¬ 
tion  of  the  pressure  doublet  in  an  arbitrary  motion  is  given.  In  order  to  sa¬ 
tisfy  the  wake  condition  it  is  assumed  that  the  pressure  doublets  move  along 
the  wake  surface  instead  of  along  the  actual  tracks  of  blades.  By  adding  the 
moving  pressure  doublets,  an  integral  equation  of  the  three-dimensional  non¬ 
steady  compressible  fluid  with  superior  singularity  is  obtained  when  the 
blades  are  in  complex  motion.  The  significant  effect  of  compressibility  on  the 
higher  harmonic  pressures  is  shown  in  this  equation.  The  kernel  function  for 
the  integral  equation  is  divided  into  two  partsi  one  with  superior  singularity 
and  a  continuous  function  with  a  strong  peak.  For  the  former  the  Hadamafd 
principal  value  can  be  determined  as  for  a  wing.  With  the  aid  of  a  proper 
spline  function  procedure  the  function  e1*0’  is  separated  from  the  function* 
hence  it  becomes  possible  to  calculate  the  higher  harmonic  pressures.  A  simple 
typical  example  is  completed  on  a  small  computer  to  justify  the  equation  and 
the  method. 
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APPLICATION  OP  KALMAN  FILTERING  TECHNIQUE  TO  AERODYNAMIC  DERIVATIVES 
FOR  A  HELICOPTER 


(Flight  Test  Research  Institute)  Yang  Songshan 

i 

ABSTRACT 

The  major  object  of  this  paper  is  to  introduce  a  method  to  obtain 
aerodynamic  derivatives  of  helicopters  using  low  pass  filtering, 

Kalman  filtering  and  least  square  techniques.  The  characteristics  of 
this  method  are  the  effect  of  high  frequency  component  of  the  rotor  is 
minimized  by  low  pass  filtering,  the  simultaneous  acquisiton  of  the 
measuring  noise  and  duration  noise  during  data  acquisition  followed  by 
Kalman  filtering  to  reduce  the  helicopter  noise  level  in  the  experiment¬ 
al  data,  and  finally  the  aerodynamic  derivatives  of  the  helicopter  are 
obtained  using  the  least  square  method.  In  order  to  increase  the 
accuracy  of  Kalman  filtering,  the  least  square  method  is  used  to  obtain 
aerodynamic  derivatives  of  the  helicopter  from  experimental  data  to  be 
used  as  the  initial  derivatives  in  Kalman  filtering.  The  calculated 
results  indicate  that  this  method  can  significantly  reduce  the  noise 
contained  in  the  experimental  data  which  cuts  ?Q%  of  the  error.  The 
computation  work  load  is  far  less  than  taht  of  maximum  similarity  method 

SYMBOLS 

Vv,V  ,V — velocity  components  along  x,y,z  axis,  respectively, 
x  y  z 

u)x,o)y,a)z — angular  velocity  components  along  x,y,z  axis,  respectively. 

ao,al,bl — coninS  angle,  reverse  angle  and  incline  angle  of  the  rotor 
respectively 

<|>^ — total  distance  of  the  rotor 
<|>wj — blade  distance  of  the  rear  rotor 

X — longitudinal  incline  angle  of  automatic  incliner 
n — transverse  incline  angle  of  automatic  incliner 
0 — angle  of  elevat-on  and  depression 
Y — angle  of  inclination 


F,  D — represent  stability  derivative  matrix  and  control  derivative 
matrix,  respectively 

w,  v — represent  duration  noise  and  measurement  noise,  respectively. 

At — sampling  interval 

E[X] — mathematical  expectation  of  X 
T 

X  — transformation  matrix  of  X 
X"1 — inverse  matrix  of  X 
x  — estimated  value  of  X 
I — unit  matrix 

received  in  May  1981 

1 .  INTRODUCTION 

Any  object  has  its  special  characteristics.  Similar  to  other 
flying  vehicles,  a  helicopter  has  its  own  characteristics  and  a  range 
of  conditions  for  best  performance.  A  helicopter  may  have  one  or  two 
rotors  and  very  high  degrees  of  freedom  in  motion.  The  experimental 
data  also  contain  higher  levels  of  noise.  Under  ordinary  conditions, 
the  longitudinal  long  period  is  basically  unstable.  These  problems 
brought  a  certain  degree  of  difficulty  in  the  processing  of  experimental 
data.  Some  of  the  methods  are  applicable  to  fixed  wing  aircraft  but 
cannot  be  used  for  helicopters.  Since  1970,  the  Kalman  filter  tech¬ 
nique  and  the  maximum  similarity  method  b&ve  been  used  to  process 
helicopter  data  In  other  countries. 

In  this  paper,  through  low  pass  filtering  and  Kalman  filtering, 
the  noise  contained  in  the  data  is  minimized.  Finally,  the  aero¬ 
dynamic  derivatives  of  the  helicopter  are  extracted  accurately  using 
the  least  square  method.  In  order  to  achieve  the  highest  degree  of 
accuracy  possible,  the  initial  derivatives  during  Kalman  filtering 
were  obtained  by  directly  applying  the  least  square  method  to  the 
original  experimental  data.  The  calculated  results  show  that  accuracy 
has  been  significantly  Increased  and  the  errors  of  main  parameters 
were  reduced  by  70%.  As  compared  to  the  method  involving  the  use  of 
least  square  method  based  on  data  obtained  from  Kalman  filtering 
reported  in  [1],  a  low  pass  filtering  loop  was  added.  The  initial 
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derivatives  of  Kalman  filtering  were  obtained  from  direct  application 
of  least  square  to  the  experimental  data  instead  of  engineering  pro¬ 
jections.  Therefore,  the  accuracy  is  higher  than  those  of  the  least 
square  method  and  the  least  square  method  using  data  based  on  Kalman 
filtering  as  described  in  [1].  This  method  requires  far  less  work 
load  in  computation  as  compared  to  the  maximum  similarity  method  and 
has  all  the  advantages  of  the  above  methods.  The  high  accuracy  in  cal¬ 
culation  and  the  simplicity  in  computation  allows  the  automatic  process¬ 
ing  of  data  as  well  as  manual  computation. 

II.  EQUATION  OF  MOTIONS  OF  THE  HELICOPTER 


ft 


ft 


ft 


In  order  to  describe  the  motion  of  a  helicopter,  at  least  six 
degrees  of  freedom  are  required.  For  a  single  rotor  helicopter,  it  § 

can  be  described  based  on  a  nine  degree  of  freedom  model. 


The  following  is  to  establish  the  nine  degree  of  freedom  equations 
for  a  single  rotor  helicopter.  For  ease  of  derivation,  let  us  assume:  ft 

(1)  the  body  of  the  helicopter  is  an  absolute  rigid  body  which  has 
six  degrees  of  freedom; 

(2)  the  rotor  of  the  helicopter  has  three  degrees  of  freedom,  aQ,  • 

al»  b1; 

(3)  the  helicopter  has  a  plane  of  symmetry,  i.e.,  longitudinal  plane. 


Using  the  body  axis  as  the  reference  coordinate  system,  the  origin 
of  the  coordinate  system  is  located  at  the  center  of  gravity.  The 
x-axis  points  forward,  y-axis  upward  and  z-axis  toward  the  right. 

Based  on  the  force  and  torque  balancing  principles,  the  following  incre¬ 
ment  function  can  oe  written 
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Rewriting  the  system  in  a  matrix  form,  using  X  to  represent  the 
change  in  status  and  U  as  the  control  input,  then  we  get 

Ax-F-AJf +f?-A£/+  w 

The  measurement  equation  can  be  expressed  as 

Z  '"H  •  &X  +  v 

If  all  the  parameters  can  be  measured  directly,  then  the  measurement 
function  is  changed  to 

Z-bX+v 

III.  NOISE  TREATMENT  AND  DERIVATIVE  IDENTIFICATION  IN  THE 
EXPERIMENTAL  DATA  OF  HELICOPTERS 

The  experimental  data  of  a  helicopter  contain  the  following  types 
of  noise: 

(1)  measurement  noise;  including  the  instrumentation  error  on 
board  and  system  error 

(2)  duration  noise;  all  the  factors  influencing  the  equations  of 
motion  will  cause  duration  noise.  For  example,  gust  disturbance, 
error  in  the  simulation  of  the  degree  of  freedom  of  the  rotor,  error 


of  nonlinear  simulation  and  system  error  on  board 
(3)  coordination  error- and  other  errors. 


In  order  to  improve  the  accuracy  of  the  experimental  results,  the 
pilot  was  requested  to  operate  precisely  to  minimize  coordination  error. 
If  possible,  the  test  was  carried  out  under  stable  airflow  conditions. 

In  order  to  eliminate  the  effect  of  the  high  frequency  portion  of  the 
rotor,  a  low  pass  filtering  method  was  first  used  on  the  experimental 
data.  Then  it  was  followed  by  the  Kalman  fitting  technique  to  minimize 
the  onboard  noise  contained  in  the  data.  Finally,  the  aerodynamic 
deviatives  were  extracted  from  the  data  using  the  least  square  method. 
The  processing  procedures  are  as  shown  in  Figure  1. 


Figure  1.  Functional  diagram 

1 — low  pass  filtering;  2 — Kalman  filtering;  3 — least  square 
method;  4 — least  square  method 


In  order  to  prevent  the  loss  of  meaningful  effective  signals  during 
low  pass  filtering,  the  cutoff  frequency  of  the  low  pass  filter  should 
be  higher  than  the  frequency  of  the  effective  signal.  In  addition, 
the  choice  of  the  cutoff  frequency  of  the  low  pass  filter  should 
satisfy  the  sampling  theory,  i.e.. 


Let  us  assume  that  Z^is  the  input  of  the  low  pass  filter  and  X. 
is  the  output,  then  the  first  order  low  pass  filtering  equation  is 

A%-(  1  - 2xf.At)Xi,.l  +  2nf.\;Zt 


Assuming  that  the  measurement  noise  and  duration  noise  are  both 
white  noise,  then  the  average  square  value  of  the  measurement  data 
should  be  the  sum  of  the  average  square  values  of  the  effective  signal 
and  the  noise  signal.  For  a  state  variable  X 


where  °> — noise  average  square  value 

— effective  signal  average  square  value 
In  addition 

o* -/./?„ 

Ra — intensity  of  noise 

If  the  same  signal  passed  through  two  filters  of  different  cutoff 
frequencies,  then  the  above  equation  can  be  written  as 

qJi-/>£-+q* 
oW.tf.  +  oJ 

thus 

Ra17zTT 


The  measurement  noise  contained  in  the  experimental  data  after 

passing  through  a  filter  with  a  cutoff  frequency  f  is 

c 


qJt-q*. , 


Similarly,  for  duration  noise,  its  intensity  is 

1  » 

a.- 

the  duration  noise  contained  in  the  data  is 

Q-Qj/.-Auft  /. 

/»—/ 1 

where 


Passing  the  filtered  data  from  a  filter  with  a  cutoff  frequency 
of  fQ  into  the  Kalman  filter,  then  the  onboard  noise  contained  in  the 
data  is  reduced  to  a  minimum. 


Since  we  assumed  that  w  are  white  noise  and  »  are  indepen¬ 
dent,  then  for  discrete  condition 
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The  Initial  state  estimate  and  the  coordination  difference 
matrices  are  x*,  P%  respectively:  _ 

i»*-£c<x*-x.)(i,-xtn 

The  observation  model  is  assumed  to  be 

Z*«/f*X*+o* 


Based  on  the  nonpartial  minimum  square  difference  estimation 
theory,  we  can  derive 
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In  order  to  prevent  divergent  phenomenon  in  the  Kalman  iteration 
operation,  the  initial  derivatives  of  the  model  are  obtained  based  on 
least  square  method  using  experimental  data. 

Because  of  passing  through  the  low  pass  filter  and  Kalman  filter, 
the  noise  contained  in  the  data  is  reduced  to  a  minimum.  Then,  the 
derivatives  are  obtained  using  the  least  square  method  which  brings 
about  more  accurate  results. 


i*i* 


Based  on  the  least  square  theory 


S(x-A’)*-  minimum 

For  a  measurement  equation  of  the  following  form 

X^4-x- 

let  us  multiply  both  sides  of  the  equal  sign  by  the  transformation 

matrix  of  X_  to  get 
m 

XlX^-Xlx- 

so 

A~{xix.yixix. 

IV.  USING  KALMAN  FILTERING  METHOD  TO  OBTAIN  THE  LONGITUDINAL 
DERIVATIVES  OF  CHARACTERISTIC  ROOT  OF  X-6  AIRCRAFT 

In  order  to  explore  the  validity  and  applicable  range  of  the 
Kalman  filtering  method,  that  method  was  used  to  process  the  experiment¬ 
al  results  of  longitudinal  stability  of  the  X-6  aircraft.  The  proce¬ 
dures  are  as  shown  in  Figure  1. 

The  experimental  results  indicate  that  this  method  can  reduce  the 
errors  of  major  parameters  by  over  70$.  Figures  2  and  3  provide  a  com¬ 
parison  between  results  obtained  using  this  method  and  the  least  square 
method.  From  the  figures,  it  is  clearly  shown  that  the  error  band  of 
the  Kalman  filtering  method  is  significantly  narrower. 

V.  CONCLUSIONS 

The  method  involving  passing  the  original  data  through  a  low  pass 
filter  and  Kalman  filter  then  followed  by  using  least  square  method  to 
obtain  the  aerodynamic  derivatives  of  helicopters  has  been  proven  to  be 
a  practical  method.  Because  a  low  pass  filter  is  added  to  the  loop 
before  the  Kalman  filter  and  because  the  original  derivatives  for  Kalman 
filtering  are  obtained  directly  from  experimental  data  using  the  least 
square  method,  the  accuracy  of  the  calculated  results  is  higher  than 
those  of  the  least  square  method  and  the  least  square  method  using 
Kalman  filtering  of  data  as  reported  in  [1].  The  computation  workload 
is  also  far  less  than  that  of  maximum  similarity  method.  This  method 
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Fig. 2  Phugoid  Motion  Darning  Comparison  of  Least  Square  Estimator 
and  Kalman  Filter  Estimator 

1— kilo  grams ;  2 — meters;  3 — Kalman  filtering;  4 — leant  square; 
5 — theoretical  estimation;  6 — 1/sec 


Fig. 3  ilf?*Demative  Comparison  of  Least  Square  Estimator  and  Kalman  Filter  Estimator 


1— ^kilograms;  2 — meters;  3 — least  square  method;  4 — Kalman  filtering 
method;  $ — 1/sec 
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APPLICATION  OF  KALMAN  FILTERING  TECHNIQUE  TO 
.  AERODYNAMIC  DERIVATIVES  FOR  A  HELICOPTER 


Yang  Songshan 
(.Flight  Test  Research  Institute ) 

Abstract 

This  paper-describes  a. -method. ior  extracting  the  aerodynamic  derivatives 
of  a  helicopter  from  flight  data  by  means  of  low  pass  filtering,  Kalman  fil¬ 
tering  and  least  square  techniques.  The  method  can  be  summaried  as  follows: 

a  )  The  high  frequency  effects  from  the  rotor  are  eliminated  by  low  pass 
filtering,  and  then  measurement  noise  and  process  noise  statistics  are 
obtained. 

b )  By  using  Kalman  filter,  the  random  noise  is  minimized  and  the  bias 
error  is  eliminated^  but  the  derivatives  are  not  identified  yet. 

c )  The  derivatives  are  identified  from  raw  data  by  least  square 
technique.  They  serve  as  initial  values  for  Kalman  filtering. 

d)  The  final  derivatives  are  extracted  from  Kalman  filtering  data  by 
least  square  technique^ 

The  method  presented  requires  considerably  less  computation  than  the  ma¬ 
ximum  likelihood  method  of  reference  Cl].  It  is  more  accurate  than  the 
least  square  technique  and  the  least  square  technique  with  Kalman  filtering 
data  of  reference  Cl].  This  method  cut  down  more  than  70  percent  of  error 
in  comperisoa  with  the  least  square  technique. 
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A  FEASIBLE  SIMPLIFIED  METHOD  FOR  FINITE  ELEMENT 
GRID  OPTIMIZATION 

I 

Beijing  Institute  of  Aeronautics  and  Astronautics,  Gong  Yaonan 

ABSTRACT  - 

One  of  the  major  considerations  in  using  the  finite  element  method 
to  carry  out  an  analysis  is  how  a  computational  result  can  be  obtained 
with  high  accuracy  and  least  computational  effort  for  a  specific  prob-  |- 

lem.  This  is  the  optimum  finite  element  discretization  problem  devel¬ 
oped  rapidly  in  the  past  10  years.  A  lot  of  optimization  criteria  and 
corresponding  procedures  were  proposed.  The  purpose  of  this  paper  is 
to  present  a  new  effective  method  for  improving  finite  element  discre-  | 

tization  in  order  to  avoid  the  difficulties  and  shortcomings  of  the  pre¬ 
sently  available  methods;  i.e.,  to  avoid  the  excessively  huge  computa¬ 
tional  time  In  the  pure  mathematical  method  and  the  difficulties  In  the 
batch  mode  method.  Its  program  is  simple  to  compile.  It  does  not  § 

require  special  software  for  computation  of  counter  lines  and  computer 
graphic  display  equipment.  A  feasible  direction  for  the  approximate 
optimization  of  discretization  can  be  obtained  from  an  error  analysis. 
Furthermore,  a  one-dimensional  search  is  used  to  obtain  an  approximate  i 

optimal  discretization.  The  numerical  example  shows  that  the  benefits 
are  apparent  by  using  the  method  described  In  this  paper. 

I .  INTRODUCTION  • 

When  using  the  finite  element  method,  similar  to  any  other  approx¬ 
imation  methods,  we  are  concerned  about  whether  it  is  possible  to  obtain 
accurate  and  convergent  approximate  solutions.  Usually,  two  aspects  • 

can  be  used  to  promote  the  accuracy  and  convergence  of  the  approximate 
solution  [2].  p-convergence  or  h-convergence .  Let  us  leave  alone  the 
fact  that  p-convergence  is  ineffective  under  many  conditions,  regard¬ 
less  of  which  of  the  two  methods  is  used  the  net  result  is  to  increase  • 

the  total  number  of  degrees  of  freedom  in  the  computational  model  which 
also  means  an  increase  in  computational  effort.  Is  it  possible  to  main¬ 
tain  the  number  of  characteristics  of  the  elents  unchanged  and  by 
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adjusting  the  grid  to  accomplish  an  optimal  numerical  solution  when 
the  above  mentioned  conditions?  Or,  in  order  to  reach  the  accuracy 
required  what  is  the  minimum  number  of  nodal  points  (number  of  degrees 
of  freedom)  increase  necessary?  Thus,  a  problem  is  laid  in  front  of 
us  which  is  the  optimum  finite  element  discretization  or  called  grid 
optimization  [1]. 

In  the  early  stage  of  study  an  optimum  finite  element  discretiza¬ 
tion,  a  mathematical  planning  method  has  been  adopted.  However,  the 
computational  effort  of  this  procedure  is  so  huge  that  the  computation¬ 
al  budget  would  not  have  been  so  high  even  if  a  uniformly  dense  grid  was 
adopted  to  obtain  the  same  degree  of  accuracy.  Therefore,  the  research¬ 
er  in  this  field  had  to  look  for  other  avenues  to  find  a  better  (not 
optimum)  finite  element  discretization  method.  In  the  recent  10  years, 
many  optimization  criteria  and  corresponding  procedures  have  been  pro¬ 
posed  [ 3—1 8 ] .  However,  these  methods,  some  of  them  have  obvious  char¬ 
acteristics  of  batch-mode  operation  and  some  of  them  require  the  graphic 
display  equipment  interfaced  with  the  computer  and  the  coordination  of 
the  special  software  which  can  compute  and  plot  counter  lines.  The 
purpose  of  this  paper  is  to  present  an  effective  method  to  improve  the 
optimum  finite  element  discretization  problem.  Under  the  condition  of 
not  increasing  the  total  number  0f  degrees  of  freedom,  the  accuracy  of 
the  solution  can  be  raised. 

II.  FEASIBLE  METHOD  SUGGESTED  IN  THIS  PAPER 

Let  the  nodal  point  coordinates  of  the  finite  element  grid  be  a 
set  of  design  variables  to  be  determined  which  is  represented  by  x.  The 
bar  under  x  indicates  that  x  is  a  vector. 

received  in  April  1981 

The  constraints  that  x  should  satisfy  are  to  assure  the  boundary 
shape  of  the  continuous  body  and  not  to  produce  irregular  elements. 

The  optimization  objective  to  be  reached  is  to  have  the  minimum  (or 
close  to  minimum)  total  potential  energy  under  an  assigned  load.  If 
the  initial  value  (nodal  point  coordinates  of  the  initial  grid)  of 
the  design  variable  x  is  x°*  The  components  to  change  x°  are 


5-«  +  fd 

where  x°  and  x  are  two  points  in  design  space.  The  vector  d  represents 
the  direction  of  motion  when  the  nodal  point  coordinates  moved  from  x° 
to  x.  t  is  the  amplitude  of  the  move.  Our  purpose  is  to  find  a  proce¬ 
dure  to  locate  a  point  x*  in  the  design  space: 

?•-?+/*$ 

which  can  make  the  total  potential  energy  to  reach  a  minimum  or  near 
minimum.  For  this  the  key  step  is  to  find  a  feasible  direction  d  using 
a  method  as  simple  as  possible.  It  is  followed  by  a  one-dimensional 
search  to  determine  the  step  length  t*. 

In  order  to  determine  the  feasible  direction  d,  we  have  to  consider 
the  error  estimation  of  piecewise  approximation  under  one-dimensional 
conditions.  Assuming  f(r)  is  the  distribution  of  a  certain  fluid  func¬ 
tion  (e.g.,  stress  energy  density)  and  the  Taylor  series  expansion  at  a 
point  r^  is 

/  (  r  )  -  /  (r,)  +  j;  *  ■ r  “r<)  V»  (r/)  +  E  ( *1 ) 

P-l 

where 

f<  71 > T,€E(r" r  > 

is  an  extra  term.  Let  us  divide  the  one-dimensional,  region  into  N-l 
subdivisions  and  each  subdivision  is  approximated  by  a  finite  element. 

In  order  words,  we  are  using  N-l  finite  elements  to  perform  piece-wise 
approximation  with  respect  to  f(r)  in  a  one-dimensional  region.  If  the 
approximate  polynomial  of  the  element  g(r)  is  of  nth  order,  then  the 
cutoff  error  of  the  ith  element  (its  two  nodal  points  coordinates  are 
^  and  r1+1)  is 

E'(n)““U+^r/<"W(T,')  n'e(r'»  r'*l) 

Therefore,  different  nodal  point  distribution  will  result  in  different 
local  error  distribution  which  in  term  causes  different  total  error. 

In  order  to  make  local  error  distributed  evenly,  we  should  let 

is, ( )-£,( *1  )  —  t|)«  constant 

where  N  is  the  total  number  of  nodal  points.  Which  means 

(r,-r.)  •V'l75rnTnJi  — —  (rw-rw.,)  '^i7T^TT0Wl  -constant 
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or 


_  constant 

'  ">'!/<-“  (n,)| 


By  adjusting  the  value  of  the  constant,  it  is  possible  to  make  r^  and 
rN  fall  on  the  two  boundaries  of  the  one-dimensional  region.  Further 
more,  equation  (1)  can  be  rewritten  as 


_  .  -  constant 


To  determine  the  nodal  point  positions  of  the  elements  from  equa¬ 
tion  (2),  we  will  obtain  a  nonuniform  local  error  distribution.  Equa¬ 
tion  (1)  is  a  special  case  of  equation  (2).  By  adjusting  the  values 
of  k  and  1,  it  is  possible  to  get  the  best  error  distribution. 


The  above  discussion  is  for  one-dimensional  conditions.  The  same 
results  can  be  applied  to  two-dimensional  situations  and  the  same  good 
results  can  be  obtained.  Although  the  corresponding  optimization  cri¬ 
teria  can  be  obtained  from  the  estimation  of  error  based  on  the  Taylor 
series  expansion  in  a  two-dimensional  situation,  yet  further  simplifica¬ 
tion  is  needed  because  an  ideal  grid  optimization  method  must  be  [8] 
matched  with  low  cost  computation  and  simple  program  design,  otherwise 
the  original  objective  of  optimization  is  not  accomplished.  The  exten¬ 
sion  of  the  conclusions  for  the  one-dimensional  case  to  a  two-dimensional 
situation  is  a  simplification  procedure.  The  actual  steps  are  in  the 
following: 


The  displacement  and  stress  energy  nodal  point  values  are  calcul 


ated  based  on  the  original  uniform  grid  *» .  With  respect  to  each  orig¬ 
inal  (both  x  and  y  directions)  grid  line,  an  extrapolated  polynomial 
f(r)  is  used  to  approximate  the  stress  energy  distribution  along  the 


grid  line: 


where  f(ri)  is  the  value  of  stress  energy  at  nodal  point  i  on  the  grid- 
line.  Then,  based  on  the  criteria  given  by  equation  (2),  the  new  coor¬ 


dinates  r^  of  the  nodal  points  (in  the  x  or  y  direction)  on  the  grid 
line  are  determined.  After  processing  all  the  grid  lines  in  both  direct¬ 


ions  using  the  method  described  above,  we  obtain  a  new  improved  grid  x 
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which  contains  the  coordinates  in  both  x  and  y  directions  of  every 
nodal  point  from 

(4) 

we  get  d’ .  The  feasible  direction  d  can  be  obtained  by  unitize  d  * . 


Sometimes  it  is  necessary  to  make  individual  corrections  for  d 
thus  created  in  order  to  avoid  the  formation  of  irregular  elements  in 
the  new  grid.  When  the  mathematical  planning  method  was  used  to  opti¬ 
mize  the  grid,  this  correction  was  contained  in- the  nonlinear  constraint 
conditions.  Since  constraints  are  not  used  here,  therefore,  corrections 
may  be  needed.  However,  in  general,  only  individual  nodal  points 
require  occasional  correction.  In  [15],  a  correction  method  was  given 
which  can  serve  as  a  good  reference . 


Whether  the  grid  has  been  optimized  can  be  evaluated  from  various 
angles  depending  upon  the  objective  of  optimization.  If  from  the  point 
of  view  of  the  entire  system  then  discretization  accomplishes  optimiza¬ 
tion.  The  potential  energy  of  the  total  system  is  a  minimum,  if  our 

% 

concern  is  focused  on  the  local  stress,  then  the  calculation  of  stress 
should  approach  its  accurate  value  as  close  as  possible.  In  order  to 
observe  the  feasibility  of  the  method  introduced  in  this  paper,  three 
examples  involving  numerical  calculations  are  given  as  below. 

III.  NUMERICAL  EXAMPLES 


As  the  first  numerical  example,  let  us  study  a  square  plate  with 
loads  concentrated  on  the  four  corners  as  shown  In  Figure  1.  Due  to 
symmetry,  only  a  quarter  of  the  plate  Is  analyzed.  Turcke  [6]  used  tri¬ 
angular  elements  and  Rosenbrock  method  to  solve  the  optimum  discretiza¬ 
tion  problem  for  this  plate.  This  paper  adopted  quadrilateral  equal- 
parametric  elements.  The  initial  even  grid  is  shown  in  Figure  1. 

Based  on  equation  (2)  and  choosing  k  *  1  ■  2,  we  obtained  the  optimized 
grid  as  shown  in  Figure  2.  The  results  are  shown  in  Table  1  where  irp 
represents  the  total  potential  energy  and  t?  indicates  the  relative 
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value  of  the  total  potential  energy  using  the  obtained  from  Turcke's 
optimized  triangular  grid  as  the  basis.  SED  represents  the  stress 
energy  density  and  NDF  represents  the  total  number  of  effective  degrees 
of  freedom. 


Turcke[6]  this  work  (v*0.33)  Carroll  (v=0.0)[l8] 


triangular  quadrilateral  f  quadrilateral 


even 


0.7601  1.3383 

40  40 

38.57  5377.46 


0.7339  1.5451 

40  40 


TABLE  1 


1.2919 

4900 


Prom  Table  1  we  can  see  that  after  the  optimization  using  the  method 
presented  in  this  paper  the  absolute  value  of  the  total  potential  energy 
can  increase  by  76$.  This  corresponds  to  an  increase  in  the  number  of 
degrees  of  freedom  from  40  to  4900.  The  benefit  is  significant. 


Fig.l  Square  plate 


litui  eren  grid 


Fig. 2  Square  plate  -  — optimized  frid 


As  the  second  numerical  example,  let  us  study  the  simple  supported 
beam  as  shown  in  Figure  3*  The  center  of  the  span  was  concentrated  with 
force.  Turcke  [17 ]  had  obtained  an  optimized  grid  based  on  the  contour 
lines  of  the  main  stress.  According  to  his  report,  the  absolute  value 
of  the  total  potential  energy  of  the  optimized  grid  increased  by  20%  as 
compared  with  that  before  optimization.  He  believed  that  this  is  a  sig¬ 
nificant  improvement.  In  this  paper,  a  four  nodal  point  isoparametric 
element  method  was  used  to  carry  out  a  calculation  with  respect  to  the 
grids  shown  in  Figures  3  and  4.  In  addition,  an  optimization  was 


TABLE  2 


grid 

even 

density 

optimized 

if 

p 

1.0 

1.0478 

1.1633 

SED  a/b 

15.19/14.93 

38.98/10.82 

3649.99/20701.91 

NDF 

43 

63 

43 

TABLE  3 


this 

method 

Carroll 

grid 

even 

optimized 

even 

optimized 

if 

p 

1.0 

1.46 

1.0 

1.46 

NDF 

12 

12 

12 

12 

a 

max 

3038.6 

4960.9 

performed  based  on  equation  (2)  by  choosing  k  =  2  and  1=0  with  respect 
to  the  grid  shown  in  Figure  3*  The  optimized  grid  is  shown  in  Figure  5. 
The  calculated  results  are  shown  in  Table  2.  In  Table  2  the  SED  has  two 
expressed  values  a/b  where  a  represents  the  stress  energy  density  at  the 
point  of  concentrated  load  and  b  represents  the  stress  energy  density  at 
the  support.  From  Table  2,  it  is  clear  that  the  method  presented  in  this 
paper  has  significant  improvement  with  regard  to  either  total  system  error 
or  partial  local  error. 

The  two  examples  above  are  on  concentrated  load.  The  third  numer¬ 
ical  example  is  a  cantilever  beam  under  evenly  distributed  load.  The 
initial  grid  is  as  shown  in  Figure  6.  According  to  equation  (2)  and 
taking  k  ■  2  and  1  ■  0,  an  optimized  grid  is  obtained  as  shown  in  Figure 
7.  The  calculated  results  are  shown  in  Table  3*  The  deflection  curve 
is  shown  in  Figure  8.  (The  numerical  computation  of  the  third  example 
was  done  by  comrade  Ah  Yejen). 

IV.  CONCLUSIONS 

The  grid  optimization  method  given  in  this  paper  is  simpler  and 
more  economical.  It  does  not  matter  whether  It  is  a  concentrated  load 
or  evenly  distributed  load  condition,  good  results  can  always  be  obtained. 
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The  computational  time,  including  the  entire  optimization  procedures, 
is  merely  5-  of  that  of  the  evenly  denser  grid  method  of  the  same 

accuracy.  This  method  was  only  explored  with  two  dimensional  balance 

problems.  Whether  it  can  be  further  extended  to  three  dimensional 
problems  is  yet  to  be  studied.  Other  problems  such  as  vibration  and 
stability  are  not  yet  dealt  with. 
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A  FEASIBLE  SIMPLIFIED  METHOD  FOR  FINITE 
ELEMENT  GRID  OPTIMIZATION 


Gong  Yaonan 

( Beijing  Insiiiuie  of  Aeronautics  and  Astronautics ) 

Abstract 

An  important  consideration  encountered  in  the  use  of  the  finite  element 
method  is  how  a  computational  result  can  be  obtained  as  accurate  as  possible 
with  the  least  computation  effort  for  a  special  problem*  especially  for  a  pro¬ 
blem  of  singularity  with  the  character  of  stress  -concentration.  This  is  the 
problem-of  Optimum-  Finite.  Element  JDi$.cretization  developed  during  the  recent 
decade.  In  this  period  a  lot  of  optimization  criteria  and  corresponding  proce¬ 
dures  were  proposed.  The  purpose  of  the  present  paper  is  to  suggest  a  new 
effective  method  for  improving  the  finite  element  discretization  which  can 
avoid  the  difficulties  and  shortcomings  encountered  in  a  batch-mode  operation 
and  doesn't  need  the  special  software  for  computation  of  contour  lines  and 
computer  graphic  displays.  After  a  relatively  simple  operation*  the  grid  opti¬ 
mization  can  be  performed  by  the  aid  of  a  feasible  direction  obtained  via  opti¬ 
mality  criteria  combined  with  the  one-dimensional  search  techniques.  The  com¬ 
putational  results  of  numerical  examples  presented  show  that  in  the  best  case 

from  an  even  grid  to  an  optimized  grid  a  gain  of  the  total  potential  energy 

*T  .  ^ 

was'  obtained^*  being  equivalent  to  -an  increase  in  the  number  of  degrees  of 
freedom  from  40  to  4900.  It  is  obvious  from  tables  and  graphics  presented 
that  quite  satisfactory  results  are  obtainable,  no  matter  whether  the  global  or 
the  local  error  is  concerned. 


THE  CRACK-FREE  LIFE  PREDICTION  FOR  STRUCTURAL  JOINTS  UNDER 
CONSTANT  AMPLITUDE  LOADS 


Aircraft  Strength  Research  Institute,  Xue  Jing  Chuan  and 
Yang  Yugong 

ABSTRACT 

This  paper  briefly  describes  the  method  to  predict  life  for  struc¬ 
tural  joints  under  constant  amplitude  loads  using  the  stress  severity 
factor  concept  and  introduced  the  procedures  to  determine  factors  a  and 
6  based  on  that  concept.  An  approach  dealing  with  the  effect  of  scratch 
is  derived  by  fitting  the  test  data.  Simultaneously,  the  reliability  of 
this  method  has  been  verified  through  the  fatigue  test  of  over  100 
riveted  joints. 

INTRODUCTION 

Structural  joints  can  be  seen  almost  in  every  part  of  the  aircraft 
structure.  How  to  determine  the  life  of  these  structural  joints  is  a 
problem  of  utmost  concern  for  aircraft  design  personnel.  This  paper 
introduces  a  stress  severity  factor  method  to  predict  the  life  of  jointed 
plate  under  axial  loads  based  on  the  S-N  curves  of  simple  notch  test 
samples  (materials). 

The  fatigue  characteristics  of  joints  are  mainly  affected  by  the 
installation  of  the  hole,  fastener  and  the  material  and  assembly  method 
of  the  matching  plates.  In  addition,  the  scratch  between  the  fastener 
and  the  hole  is  also  influencing  its  life.  The  stress  severity  factor 
Is  an  indicator  of  fatigue  quality  which  takes  the  above  factors  into 
consideration.  The  effect  of  scratch  is  considered  through  the  correct¬ 
ion  from  the  S-N  curve  of  the  material.  It  should  be  pointed  out  that 
this  paper  is  limiteded  to  a  one-dimensional  stress  problem. 

This  paper  has  been  reviewed  by  Chief  Engineer  Fong  Chungyuch.  The 
work  was  assisted  by  Associate  Professor  Yang  Chingsbon  of  Northwest 
Industrial  University.  Assistance  also  came  from  comrades  of  Shenyang 


Aix „^aft  Company  and  Hongan  Aircraft  Company.  Sincere  appreciation  to 
those  who  assisted  is  expressed  here. 

I.  STRESS  SEVERITY  FACTOR  [1] 


Stress  severity  factor  is  a  parameter  which  reflects  the  effect 
of  the  installation  of  the  hole  on  the  fastened  plate  as  well  as  the 
materials  and  the  assembly  method  of  the  fastener  and  the  matching  plate. 


As  shown  in  Figure  1,  the  maximum  local  stress  along  the  edge  of 
the  hole  is 

a—~K"-§r*  +K'-wi 

the  combined  stress  concentration  factor  is 


In  order  to  indicate  the  effect  of  surface  status  and  assembly  type 
under  fatigue  load  conditions,  two  coefficients  a  and  8  are  introduced 


to  the  define  stress  severity  factor  as: 


received  in  December  1980. 

Kr~  ~wi)  ( 1 ) 

where  R — transfer  load  of  the  fastener 

— local  side  load 
P 

— theoretical  stress  concentration  coefficient 
Kjy — compressive  stress  concentration  coefficient 
a — hole  surface  coefficient 

8 —  hole  filling  coefficient 

9 —  compression  distribution  coefficient 

In  the  elastic  region,  equation  (1)  can  be  rewritten  as: 

9  )+/,(*, )5 

Therefore,  Ky  is  only  affected  by  the  structural  form  (a,8,9,Kjy,Kt ) 
and  not  influenced  by  the  properties  of  the  material  and  stress  load. 
When  the  material  and  fastener  enter  the  modern  region,  this  is  no 
longer  true. 
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Fig.l  Stress  distribution  at  the  fastener  hole 


Fig .  2  «  and  P  specimen s  - 

i — plate;  2 — washer;  3 — screw  head;  4 — crack  warning  transducer 

The  KyO,t  at  stress  concentration  point  in  the  structure  approx¬ 
imately  represents  its  life.  Therefore,  it  is  possible  to  obtain  the 
desired  uniform  structural  life  by  the  proper  design  so  that  this  value 
is  as  identical  and  as  low  as  possible.  In  the  prediction  of  structural 
life,  K  corresponds  to  the  theoretical  stress  concentration  coefficient 

V 

II.  a,  8  COEFFICIENTS 

This  paper  does  not  wish  to  discuss  the  load  distribution  compu¬ 
tation  of  Joints.  Therefore,  the  remaining  major  problem  in  the  deter¬ 
mination  of  stress  severity  factor  is  to  determine  a  and  8.  On  the 
basis  of  the  stress  severity  factor  concept,  this  paper  proposes  the 
following  method  to  determine  these  coefficients: 


1.  Determination  of  a  coefficient 


Two  samples  with  boles  made  of  two  different  materials  and  identi¬ 
cal  geometric  shape  were  used  as  the  sample  (see  Figure  2).  The  theor¬ 
etical  stress  concentration  coefficient  is  Kt .  One  sample  has  a  stand¬ 
ard  hinge  hole,  the  other  has  a  drilled  bole.  Let  us  assume  for  the 
hinge  hole  o  ■  1.  Based  on  equation  (1),  because  3*  1,/?—  o,  Pp/Wi^a^ 
therefore, 

K„-K, 

where  is  the  stress  severity  factor  under  the  hinge  hole  condition. 
Similarly,  the  stress  severity  factor  the  drilled  hole  is 

Kr2-aK, 

Using  a  certain  level  of  stress  for  the  fatigue  test  of  the  drilled 

hole,  its  life  is  N.  .  On  the  S-N  curves  of  the  standard  hinge  hole 

with  various  stress  concentration  coefficients,  it  is  possible  to  find 

the  theoretical  stress  concentration  coefficient  Ktl  which  corresponds 

to  the  same  level  of  stress  and  life  as  N  : 

a 

Krz~K,t~aK,  (2) 

a  =Ktl/K. 

2.  Determination  of  the  B  coefficient 


The  geometric  shape  of  the  specimen  is  identical  to  that  with  a 
hole.  The  only  difference  is  that  inside  the  hole  there  is  a  non-load 
conducting  fastener.  According  to  equation  (1),  the  stress  severity 
factor  of  a  screwed  specimen  is 

Using  a  certain  stress  level  to  conduct  the  fatigue  test  of  a 
screwed  specimen,  its  life  is  Ng.  The  theoretical  stress  concentration 
coefficient  K,t  under  similar  conditions  in  terms  of  stress  level  and 
life  as  Ng  can  be  found  from  the  S-N  curves  of  standard  hinge  hole  of 

various  stress  concentration  coefficients.  Then 

Kru-a8K,-K,t 

The  results  are  shown  in  Table  1. 


(3) 


Table  1 


: 

i  i 

l  a 

1  9 

/  j 

tt  ft  1 

*  Sf  x 

1 

i 

\gf  XKCO 

l 

^  tt  ft 

mum 

1.07*  j  | 

7  xitco  j 

!  11  i 

3.  . 

j  !  1*  *** 

0.86 

|  f  3t*CO 

0.75—0.9 

1 — hinge  hole;  2 — drilled  hole;  3 — screw  connection;  4 — this  study; 

5 — reference  [4];  6 — this  study;  7 — reference  [4];  8 — this  study; 

9 — reference  [4] 

III.  THE  EFFECT  OF  SCRATCH  [2] 

Under  fatigue  load  situation  the  fastener  under  load  moves  relative 
to  the  assembly  hole  to  scratch  the  hole  wall  which  reduces  the  fatigue 
life  of  the  basic  plate.  Based  on  the  experimental  data  of  joints,  we 
gradually  obtained  the  relation  between  the  S-N  curves  with  correction 
for  the  affect  of  scratch  and  the  S-N  curves  of  materials  with  a  notch 
by  fitting  the  data: 


high  load 
mediun  load 

lg^v„.  =.  0 . 9446lgiV^. -r  o .  1109 

1 

(4) 

low  load 

0. 7686  lg<V^  +  1.0217 

S 

where  N,»  T-life  after  correction  for  scratch 
— life  of  material  with  a  notch 

The  above  relation  was  obtained  from  the  S _ -N  curve  corresponding 

max 

to  a  stress  ratio  R.  Based  on  this,  it  is  possible  to  give  the  con¬ 
stant  life  curve  of  the  joint  with  correction  for  the  effect  of  scratch 
(see  Figure  3). 
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10  30 

S.(£Jr/**K)  3 


Fig. 3  Constant  life  eurre  taking  account  of  cratch 

2  2 
1 — (kg/mm  );  2 — material;  3 — (kg/mm  ) 

IV.  PROCEDURES  TO  DETERMINE  THE  LIFE  OF  JOINTS  UNDER  CONSTANT 
AMPLITUDE  LOAD  USING  THE  STRESS  SEVERITY  FACTOR 

1)  Determine  the  internal  force  distribution  of  the  joint  under ^ 

fatigue  load  in  order  to  determine  the  transfer  load  of  the  fastener, 

A 

the  side  load  P  and  the  reference  stress  o.*  ,  etc. 

Sr 

2)  Base  on  the  geometric  shape  and  assembling  form  the  joint  to 
determine  3,  AT*,  AT,,  0,  that  the  stress  severity  factor  K  at  all 

*r 

the  stress  concentrating  points  are  known. 

3)  Determine  the  position  of  maximum  product  of  stress  severity 
factor  and  reference  stress  Kra^ , 

4)  Base  on  the  stress  severity  factor  and  reference  stress  at  that 
point  to  locate  the  predicted  life  by  looking  into  the  corresponding 
life  curves  of  identical  material  and  equal  stress  concentrating  coeffi¬ 
cient  with  the  affect  of  scratch  taken  into  account. 


61 


iam 


**LYI2CZ  #»«JLY12CZ  ^*TLY10-3Y- «-  M 


201  201  20  .12)12  20  1  20  20 


«SlY12CZ 

-L _ 1 


- -  _ 

«aHSLYl2CZ 


graEi 


Fig.S  Specimen  2 


raiaal 


lanii 


Fif  8  S«t  aftiosr  beading 

1,2— crack  point;  3 — base  plate;  4 — joint  plate;  5 — rivet;  6 — rivet  heads  on  same 
side  of  base  plate;  7,8 — typical;  9,10— crack  point;  11 — typical;  12 — constant  space 
array;  13— base  plate;  14— joint  plate;  15 — rivet;  16 — rivet  beads  all  on  same  side 
of  base  plate;  17— organic  glass;  18,20— oil  paper;  19 — specimen;  21 — wooden  plate 
22— plastic  foam. 
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V.  JOINT  FATIGUE  TEST 


1)  The  specimens  have  two  forms,  i.e.,  the  rivets  are  aligned 
parallelly  and  mismatched  .  The  basic  plate  thickness  is  4  mm.  The 
joint  plate  thickness  is  5  mm.  The  plate  material  is  LY12-CZ.  The 
rivets  nave  burled  heads.  The  diameter  is  6  and  material  is  LY10. 

The  fiber  direction  of  the  plate  material  is  parallel  to  that  of  stress 
The  shapes  of  the  specimen  are  shown  in  Figures  4  and  5. 


2)  The  specimens  are  sandwiched.  In  order  to  effectively  reduce 
the  side  vibration  and  bending  of  the  specimen  under  load,  a  guidance 
plate  was  installed  as  shown  in  Figure  6. 


3)  The  static  destruction  was  shearing  of  rivets.  Fatigue  des¬ 
truction  all  occurred  at  the  first  row  of  rivets  on  both  sides.  The 
experimental  results  are  shown  in  the  S-N  curves  in  Figures  7-11 • 

VI.  PREDICTION  OF  JOINT  LIFE  [3,4] 

1.  Calculation  parameters 


Let  us  choose  a 
of  rivet  head). 


1.1;  0  »  0.75,  9  *  1.4  (considering  the  effect 


Kjy  *  1.3  for  parallel  array 

■  1.3  for  mismatched  array  with  four  boles  in  a  row 

*  1.2  for  mismatched  array  with  three  holes 

for  parallel  array  *  3-02  for  the  center  hole 

■  3.2  for  the  holes  on  the  side 
for  mismatched  array  *  3.2  for  4  boles  in  a  row 

*  3-1  for  3  boles  in  a  row 


2.  Conversion  between  K  and  S 

O' 

In  the  calculation  of  stress  severity  factor,  corresponding  to  a 
fatigue  load  level  (Pmax  and  P  ^  ),  the  stress  severity  factors  are 
Kymax  and  Kymln  respectively.  In  the  elastic  region,  Kymax  -  Kymln. 
If  the  fastener  enters  a  molding  stage,  then  the  two  are  not  equal. 


2  3  4  S  t  1 


Fig. XI  S-N  cam 


1 — experimental;  2 — calculated;  3 — parallel  array;  4 — (kg/min  ) 


In  order  to  carry  out  the  calculation  for  life  prediction,  the  stress 

severity  factor  Ky  *  Kymax  corresP°ndinS  t0  Pmax  and  Pmin  is  u5ed* 
Then  the  equivalent  stress  level  is: 

Kv~KymaM  | 

5-- (5^+51.)/  2  > 

2  ) 


where 


S'm  j.  ”  Ky 51  In/ Kfm\n 


Then,  life  can  be  calculated  based  on  the  procedures  described  in 
section  IV  of  this  paper.  The  results  are  shown  in  Figures  7-11* 

VII.  DISCUSSION  OF  RESULTS 


1)  The  stress  severity  factor  method  was  orignally  proposed  by 
Jarfall  [1].  It  was  further  developed  and  used  in  some  countries  in 
the  fatigue  design  of  aeronautic  structures.  This  study  considered  the 
modeling  effect  of  the  rivet.  In  this  paper,  the  stress  severity  factor 


concept  was  used  as  the  basis  to  obtain  coefficients  a  and  6  and  the 
results  were  also  given.  This  paper  also  proposed  an  empirical  method 
to  determine  the  effect  of  scratch  through  experiments.  This  paper 
extended  the  stress  severity  factor  method  to  the  calculation  of  the 
S-N  curve  of  joints.  It  provided  a  base  for  the  preliminary  life  pre¬ 
diction  calculation  under  programmed  load  conditions. 

2)  In  order  to  verify  the  reliability  of  this  method,  a  batch  of 
riveted  joints  were  fatigue  tested  and  the  following  conclusions  were 
obtained : 

The  calculated  destructed  portion  coincides  with  the  experiment. 
The  ratios  between  calculated  results  and  experimental  results  are 
mostly  between  1.0  and  0.75* 

This  method  was  further  verified  in  the  fatigue  experiment  and 
results  analysis  of  a  certain  wing  joint  of  an  aircraft. 

3)  There  must  be  a  lot  of  shortcomings  and  errors  in  the  problems 
described  in  this  paper.  Your  comments  are  welcome. 
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THE  CRACK-FREE  LIFE  PREDICTION  FOR  STRUCTURAL 
JOINTS  UNDER  CONSTANT  AMPLITUDE  LOADS 

Xue  Jingchuan  and  Yang  Yugong 
C Aircraft  Sirengih  Research  Insiiiuie') 

Abstract 

This  paper  briefly  describes  a  method  of  detailed  designing  and  predicting 
the  fatigue  life  for  structural  joints  under  constant  amplitude  loads  by  means 
of  the  Stress  Severity  Factor  concept 

J- 8 +*'  ii) 

Based  on  the  Stress -Severity  Factor  concept*  an  approach  for  determining 
factors  °  and  3  is  developed  and  some  of  their  values  are  derived  from  fa¬ 
tigue  tests. 

An  approach  dealing  with  the  cratch  is  created  by  fitting  test  data,  and 
the  relation  between  the  corrective  S~N  curves  of  structural  joints  and  the 
S~N  curves  of  materials  is  defined. 

More  than  100  specimens  of  revited  joints  were  tested.  The  test  data 
obtained  were  compared  with  the  analytical  results  favorably  and  thereby 
testified  the  feasibility  of  this  method. 

In  addition,  the  present  method  was  employed  to  predict  the  full-scale 
S-N  curves  of  structural  joints  for  their  .preliminary  design  and  gave  quite 
good  results  under  the  given  conditions. 
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AN  ANALOGY  METHOD  FOR  CRACK  INITIATION  LIFE  PREDICTION 


Zhang  Fuze 

Chinese  People’s  Liberation  Armed  Forces,  Air  Force  Research  Institute 
ABSTRACT 

This  paper  presented  an  analogy  method  for  crack  initiation  life 
prediction  of  components.  It  employs  the  life  (or  useful  life)  of  a 
known  component  obtained  from  endurance  tests  under  spectrum  load  to 
predict  the  life  of  some  type  of  components  untested.  The  prediction 
formula  of  this  analogy  method  was  deduced  in  this  paper.  It  could 
practically  eliminate  the  scattered  effect  of  the  constant  Q  in  the  Miner 
equation  ^2-^—0^  •  Therefore,  it  was  capable  of  improving  the  accur¬ 
acy  of  the  life  prediction.  This  paper  also  gave  the  analogy  prediction 
of  specimens  under  five  types  of  spectrum  load  and  an  example  of  analogy 
prediction  of  large  scale  structure  components  of  an  airplane.  It  was 
also  verified  by  fatigue  tests. 

In  the  research  of  analogy  life  prediction  in  1979,  Professor  Kao 
Chengtung  of  Beijing  Aeronautics  and  Astronautics  Institute  has  given 
his  guidance.  The  verification  work  was  carried  out.  by  Engineer  Ho 
Tungen  of  the  Air  Force  Research  Institute.  In  the  analogy  computation 
of  specimens.  Engineers  Koo  Ming  yuan  and  Zhang  Shecbak  provided  exper¬ 
imental  data  without  reservation.  Many  thanks  to  these  people. 

I .  INTRODUCTION 

As  Is  commonly  known,  ever  since  Miner’s  theory  was  established 
in  19**5,  although  it  is  still  widely  used.  It  still  has  some  inadequacy 
and  inaccuracy.  The  most  serious  problem  Is  the  scattering  of  the  con¬ 
stant  Q  in  the  equation.  In  theory  Q  -  1J  however,  in  practice  it 
scattered  quite  widely  (several  times  in  several  10  times  in  difference). 
This  brings  significant  errors  in  life  prediction  calculation. 

This  paper  deduces  an  analogy  crack  initiation  life  prediction 
calculation  equation  of  components  based  on  the  Miner  theory  and  the 
equation  SmN—c  .  In  the  computation,  it  is  not  necessary  to  manually 
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select  the  highly  scattered  Q  value.  Therefore,  it  is  possible  to 
raise  the  accuracy  and  reliability  of  life  predict!'  Consequently, 
full  scale  fatigue  tests  can  also  be  reduced. 

In  the  design  of  airplanes  and  the  calculation  of  airplane  char¬ 
acteristics,  the  analogy  method  has  been  used  over  a  long  period  of  time 
and  is  widely  applied.  However,  using  an  analogy  method  to  calculate 
the  fatigue  life  of  components  has  not  reached  a  mature  stage  with  methods 
and  complete  set  of  equations.  The  author  of  this  paper  derived  an 
earning-loss  analogy  equation  (same  principle  as  the  one  used  in  this 
paper)  in  1969  and  carried  out  some  practical  applications. 

II.  THE  ESTABLISHMENT  OF  ANALOGY  LIFE  PREDICTION  EQUATION 

FOR  CRACK  INITIATION 


Use  the  Miner  equation  to  start  with: 

K 


(1) 


where  \ — cycle  number  of  components  with  known  life 

ni — the  cycle  number  of  the  ith  order  load  in  one  cycle  of 
spectrum  load  for  a  component  with  known  life 


N^ — destruction  cycle  number  under  ith  order  load  and  then  /52 

introduce  the  equation 

S7N-c  (2) 

where  S& — stress  amplitude  of  the  load 

N — destruction  cycle  number  of  the  material  corresponding  to 
a  stress  amplitude  of  Sa. 

m  and  c — relevant  constant  characteristics  to  the  material, 
shape  and  cyclic  load  of  the  component. 

From  equation  (2),  for  the  same  component  at  the  same  location, 
there  is 

SZIN'-S:^,  (3) 
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i.e . , 

N,  _(S 

m;  Kst) 

where  Sap — an  arbitrary  constant  stress  amplitude 

Np — destructive  cycle  number  corresponding  to  Sap. 

Now,  let  us  multiply  equation  (1)  by  N  , 

K  P 


Prom  equations  (4)  and  (5): 


X  v  n*U.-QN9 
iTi  Nl 


Similarly,  the  expression  for  the  components  with  unknown  life  is 

i/  V  S'  W  (7) 


Under  the  conditons  that  the  structural  form,  load  properties, 
loading  sequence  and  material  type  are  basically  the  same  for  components 
with  known  and  unknown  lives,  we  have 


Ms 


Equation  (9)  is  the  computational  analogy  equation  for  the  predic¬ 
tion  of  crack  initiation  life  of  components.  If  the  spectrum  load  is 
not  a  stress  spectrum  but  an  overload  spettrum,  then  it  is  expressed 
in  the  following. 

In  the  elastic  region,  for  the  same  component  at  the  same  location, 
the  stress  amplitude  is  proportional  to  the  overload  increment.  Then 
we  have 
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X' 


X 


(10) 


&<tr 


where  Ag^ — the  1th  order  overload  Increment  of  component  with  known  life 
Ag.* — the  itb  order  overload  increment  of  component  with  unknown 
life 

Agp — arbitrary  overload  increment 

Equation  (10)  is  the  analogy  calculation  equation  for  known  over¬ 
load  spectrum.  In  equations  (9)  and  (10),  because  Sap  and  Agp  are  arbi¬ 


trary  quantities,  then  let  S, 
tions  (9)  and  (10). 


1  or  Agp  *  1  to  further  simplify  equa- 


As  for  the  analogy  prediction  for  identical  components,  in  engineer¬ 
ing  it  is  possible  to  use  the  average  value  m  as  an  approximation  and 
further  to  eliminate  the  effect  of  cyclic  characterics  on  the  value  of 
m  to  make  ra^  *  *  m  and  to  convert  the  stress  amplitude  to  stress 

increment  (AS  »  Simax_Simin^ *  Thus,  the  following  simplified  equation 
can  be  obtained: 


K 

S»/AS7 

-T1 - * 

2<Asr 

»-i 


jc 

X'--U! - x 

<-1 


(id 


(12) 


Based  on  the  above  equations,  as  long  as  the  value  of  m  is  obtained, 
the  unknown  life  X’  can  be  determined. 

III.  VALUE  OP  ffl  AND  ITS  EFFECT  ON  LIFE 
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1.  The  value  of  m  during  the  crack  initiation  period  is  usually 

determined  experimentally  in  engineering  where  good  s-N  curves  are  not 

available.  Under  the  conditions  that  no  test  was  performed,  it  is 

usually  taken  as  5-10  [1].  Prom  the  analogy  equation,  the  error  of  the 

relative  life  obtained  from  the  analogy  method  using  m  is  smaller  than 

that  of  the  absolute  life  obtained  directly  from  the  Miner  equation  by 

choosing  Q.  This  is  because  whether  the  m  value  was  too  large  or  too 

small,  the  numerator  and  the  denominator  would  increase  or  decrease 

1 

simultaneously.  The  effect  on  the  relative  life  (X  /X)  obtained  from 
the  analogy  method  is  very  small.  For  example,  the  data  and  curve 
shown  in  Figure  1  were  obtained  based  on  different  m  values  using  the 
analogy  method.  They  were  calculated  based  on  the  identical  components 
of  two  different  airplanes  under  their  own  spectrum  load  using  the 
analogy  method. 


Fij.  I  Relation  between  A'/X  with  n» 


From  Figure  1  we  can  see  that  the  relative  life  decreases  with 
increasing  m  value,  but  the  change  is  gradual.  When  the  value  of  m  is 
increased  by  1,  the  relative  life  is  decreased  by  51.  In  the  Miner 
equation,  when  the  constant  Q  is  increased  by  1,  then  the  life  is  in¬ 
creased  by  100$.  Furthermore,  the  variation  of  Q  in  the  calculation 
is  large  which  causes  a  larger  error  in  the  calculation.  Based  on 
these  results,  we  know  that  for  the  sample  component,  the  analogy  method 
will  not  cause  large  error  in  life  prediction  due  to  er  *or  in  choosing 
m.  Therefore,  this  analogy  method  not  only  eliminated  the  effect  of  Q 
value  on  the  Miner  equation  but  also  demonstrated  its  superiority.  In 
addition,  its  superiority  has  been  demonstrated  by  the  reduction  of  life 
prediction  due  to  the  decrease  in  the  error  caused  by  the  m  value. 


2.  The  m  value  during  the  crack  initiation  stage,  when  perfect 

» 

S-N  curves  are  availble,  can  be  used  to  obtain  and  und'er  differ¬ 
ent  cycling  characteristics  to  further  improve  the  accuracy  of  the 
estimation.  Under  this  condition,  it  is  also  possible  to  directly  use 
Miner's  equation  and  eliminate  Q  to  perform  the  analogy  procedures. 

IV.  ACTUAL  SAMPLE  OP  ANALOGY  METHOD  FOR  PREDICTION  OF' CRACK 
INITIATION  LIFE 

1.  Specimen  analogy  method  example 


The  experimental  fatigue  data  used  in  the  specimen  forming  life 
prediction  using  the  analogy  method  were  obtained  by  Beijing  Materials 
Research  Institute  in  1980.  The  loading  spectra  used  in  the  tests  were 
variable  average  spectrum,  constant  average  spectrum,  double  wave  (two- 
flow)  spectrum,  fatigue  indicator  spectrum  and  constant  damage  spectrum. 
The  tester  was  a  pcl60N  hydraulic  servo  fatique  tester  of  the  Schenck  Company. 

The  loading  sequence  of  all  the  spectra  is  low-high-low.  The  specimen 
material  is  lead  alloy  Lyl2-CZ.  Its  mechanical  properties  are  shown  in 
Table  1.  The  thickness  of  specimen  plate  is  2.5  mm  and  Kfc  =  4. 


TABLE  1.  Mechanical  properties  of  the  specimens 


o  (kg/mm2) 


(kg/mm2) 


«10(,) 


35.8 


47.8 


17.4 


The  analogy  method  was  used  for  this  specimen  using  the  fatigue 
test  results  which  give  more  accurate  programmed  spectra.  The  procedures 
of  the  analogy  method  involve  the  use  of  fatigue  test  life  of  a  set  of 
specimens  under  a  particular  spectrum  load  to  predict  the  life  of 
another  set  of  specimens  under  another  spectrum  load.  Then  the  cal¬ 
culated  results  are  compared  with  those  obtained  from  experiments  to 
determine  the  relative  error  between  the  two  conditions.  Let  us  now  use 
the  constant  average  spectrum  and  fatigue  spectrum  as  an  example  to 
carry  out  the  analogy  prediction  with  results  shown  in  Tables  2,  3  and  4. 
The  m  value  in  the  calculation  took  into  consideration  two  conditions: 
one  condition  was  to  consider  the  effect  of  the  cyling  characteristics 


R^  on  m  which  involved  the  computation  of  the  values  of  at  different 
and  the  other  condition  was  not  to  consider  the  effect  of  R^  on  m 
which  used  the  average  value  of  m  in  the  computation. 

Similar  methods  were  used  to  carry  out  analogy  calculations  to  pre¬ 
dict  specimen  life  under  variable  average  spectrum,  double  wave  spec¬ 
trum  and  constant  change  spectrum  loads.  In  the  meantime,  the  Miner 
equation  was  also  used  to  calculate  the  life  of  the  specimen  and  the 
results  were  compared  with  those  of  the  fatigue  test.  They  are  also 
shown  in  Tables  5  and  6.  In  order  to  compare  the  calculated  results 
with  the  experimental  ones  more  directly  and  realistically,  no  correct¬ 
ion  was  made  for  the  calculated  and  tested  lives.  Therefore,  the  cal¬ 
culated  lives  in  Table  5  and  6  are  true  reflections  of  the  analogy 
method  with  respect  to  the  Miner  equation.  The  fatigue  tested  life  is 
the  realistic  view  of  the  specimen  under  spectrum  load.  If  the  tested 
life  of  the  specimen  has  already  taken  scatter  coefficient  and  other 
data  treatment  into  consideration,  the  calculated  life  based  on  the  analogy 
method  contained  the  same  scatter  coefficient  and  data  treatment. 

The  S-N  curves  used  in  the  above  calculation  were  also  determined 
by  Beijing  Materials  Research  Institute.  The  properties  of  the  material 
of  the  specimen  and  Kt  are  identical  to  those  for  the  specimens  used 
in  the  above  mentioned  five  fatigue  tests.  Therefore,  the  S-N  curves 
used  in  this  analogy  calculation  are  relatively  more  accurate  and  relevant. 


2.  Analogy  prediction  example  of  fatigue  Initiation  life  of 
large  scale  structural  pieces 

The  fatigue  test  data  used  in  this  analogy  prediction  were  the 
results  of  full  scale  fatigue  experiment  of  the  main  stress  parts  on 
the  wing  and  fuselage  of  airplanes  of  the  same  type  for  different  appli¬ 
cation®.  This  fatigue  test  of  large  scale  structural  components  was 
performed  on  a  special  equipment.  The  error  of  load  increment  and  load 
indicator  is  less  than  2%. 


Table  3  The  analofy  computation  of  the  specimens 


Key:  1 — parameters  and  amputation  of  the  fatigue  indicating  spectrun;  2— average;  3 — considering  the  effect 
of;  <■  -not  considering  the  effect  of;  5— data  in  the  table  are  the  full  cycle  values  of  the  fatigue 
indicating  spectrum;  &  -not  considering  the  effect  of  R,  ,i.e. ,  m. -const  (average  m  value)  in  the 
amputation;  7— a  cycle  is  50  flight  hrs;  8— notes.  1 


Table  4  The  analogy  computation  of  the  specimens 


-«*«►<*+)  -35.0,  WT'-SS.Ox  50- 1750*  8*  'j. 

(2)-fi*»if-j*ft«T,  -«*#  <**)  »*i**»****X -64.0,  KT -84.  Ox  SO -3200*1*  r^~, 


to  ««**»****  a  -«4.o)  *tttt**j$*T«*#**x'.  **2.  * 

3»tt***,  aifctt*.  ~~  ,y-'' 


1  - 1 _ 

29 

2  "<  AflJ"** 


X  *  l®T'  -36.75x50- X837*W<MSX>?£\ 


93738.85 


(2)  ««»«»«***  (X'  - 35.0)  £ttttX«3H-»TM*#;**X,  **2.  \J_ 

*3  »«■*«*.  *£«•»,  ^ 


»-l _ 

23 

2  mtaV 


X'-  93~lla2lirJ'"62  0^  ®  "62.0  x  50  -3XOO*l*(itft  3-JO-  - 


1 — fatigue  tested  life;  2 — analogy  confuted  life  and  error;  3 — (1 )  under  con¬ 
stant  average  spectrum,  average  fatigue  test  life  A' *35.0  for  a  set  of  speci¬ 
mens  (six),  i.e.,  T '  =35 . 0x50=1 750  hrs;  4 — (2)  under  fatigue  indicating  load, 
average  fatigue  test  life  \=64.0  for  a  set  of  specimens  (six),  i.e.,  T =6 4. 0x50= 
3200  hrs;  5 — (1)  under  fatigue  indicating  spectran  load,  using  the  specimen 
fatigue  test  life  (A*64.0)  to  compute  the  specimen  life  A'  under  constant 
average  spectrun  load  vising  analogy  confutation.  Taking  the  results  in  Tables 
2  and  3  to  perform  analogy  confutation ;  6 — i.e.,  T '=36.75x50=1837  hrs  (error 
5%) ;  7 — (2)  under  the  constant  average  spectrum  load,  using  the  specimen 
fatigue  life  (A' =35.0)  to  confute  specimen  life  A'  under  fatigue  indicating 
spectrun.  Taking  the  calculated  results  in  Tables  2  and  3  to  carry  out  analogy 
_ confutation;  8 — i.e.,  T  =62. 0x50=3100  hrs  (error  3%) 


TableS  The  comparison  between  the  results  of  fatigue  tests  and  analogy 


computation 


Table  6  Specimen  lifescalculated  with  Miner's  formula 
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fatigue  indicating;  1 2 — double  wave  (2  flow);  13 — equal  damage;  14 — note? 

15 — only  the  accurate  solutions  (considering  the  effect  of  R. )  are  listed  for 
analogy  oanputation  1 


Table  7  The  comparison  between  the  lives  estimated  by  analogy  and  the 
fatigue  test  lives  of  large  components 
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Key:  1— specimen;  2— results;  3—item;  4— relative  life  ( X/X 1 )  of  the  two  structural 
components  based  ext  analogy  computation;  5 — relative  life  (X/X * )  of  the  two 
components  in  fatigue  test;  6 — relative  error;  7 — equation;  8 — overloading 
spectrum;  9 — stress  spectrun 

The  spectrum  load  in  the  fatigue  test  was  measured  using  the  same 
instrument  on  two  types  of  airplanes.  The  compilation  method  of  the  two 
spectra  Is  identical.  The  load  increasing  sequence  is  also  the  same 
(low-high  type).  Under  these  two  spectrum  loads,  the  fatigue  crack 
positions  of  the  two  structural  pieces  are  also  the  same. 


78 


In  order  to  compare  the  effects  of  overloading  spectrum  and  stress 
spectrum  on  the  prediction  equation,  the  analogy  computation  of  this 
large  scale  structural  component  was  carried  out  using  these  two  spectra. 
The  stress  spectrum  was  obtained  by  directly  measuring  at  the  cracked 

f 

portion.  The  and  m^  used  in  the  computation  are  obtained  from  the 
S-N  curves  together  with  the  consideration  of  the  effect  of  R.^.  However, 
since  the  variation  is  not  too  large,  the  average  value  m  was  used  in 
the  computation.  The  results  are  shown  in  Table  7. 

V.  CONCLUSIONS 

1.  This  analogy  life  prediction  method,  through  the  verification 
of  analogy  computation  of  the  five  sets  of  specimens  under  five  differ¬ 
ent  types  of  spectrum  load  and  large  scale  structural  component  analogy 
computation  and  fatigue  tests  has  demonstrated  that  for  the  prediction 
of  life  of  components  of  similar  type,  as  long  as  the  determination  of 
m  is  reasonable,  the  analogy  computation  presented  in  this  paper  pro¬ 
vides  results  of  higher  accuracy.  They  also  agree  very  well  with  fatigue 
test  results.  Under  this  condition,  the  effect  of  the  cycling  character¬ 
istics  on  m  (average  value  m)  can  be  neglected  and  equations  (11)  and 
(12)  can  be  used  directly. 

2.  This  analogy  computation  method  to  predict  life  is  more  con¬ 
venient  for  the  life  prediction  in  the  design  stage.  Using  the  life 
(known  life)  of  the  component  of  the  original  airplane,  the  life  of  the 
newly  designed  similar  components  can  be  calculated  using  this  analogy 
method.  Thus,  it  is  possible  to  figure  out  whether  the  newly  designed 
component  has  a  higher  or  lower  life  of  the  original  standard  airplane. 

It  is  then  very  meaningful  in  design  changes.  It  is  also  possible  to 
calculate  which  components  have  a  weak  life  and  which  components  have  a 
conservative  life.  This  serves  as  a  reference  to  the  evaluation  of  the 
life  of  the  entire  airplane. 

3.  When  the  spectrum  load  of  individual  airplanes  in  use  or  the 
spectrum  load  of  airplanes  (or  spectrum  load  of  fatigue  test)  varies, 
this  analogy  method  can  very  conveniently  compute  the  variation  of  life 
caused  by  these  changes. 
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4.  From  the  computed  results  in  Tables  5  and  6,  it  is  clear  that 
the  results  obtained  using  this  method  have  a  large  difference  as  com¬ 
pared  to  those  with  Miner's  equation.  The  causes  are  (1)  the  effect  of 
scatter  of  Q  is  reduced,  (2)  this  analogy  method  is  a  combination  of 
the  Miner  equation  and  expression  of  S-N  curves.  It  is  not  necessary 

to  repeatedly  check  the  S-N  curves  or  constant  life  lines  in  computations. 
Therefore,  the  errors  of  the  checking  and  plotting  figures  can  be  drasti¬ 
cally  reduced;  (3)  the  effect  of  m  on  relative  life  is  reduced. 

5.  For  the  analogy  computations  of  different  components  under  the 
conditions  that  the  load  characteristics  and  working  sequences  are  basi6- 
ally  the  same,  this  method  is  still  applicable  based  on  the  Miner  theory. 
However,  due  to  the  difference  in  structural  type  and  material,  a  cer¬ 
tain  degree  of  error  is  created.  Therefore,  some  correction  is  necess¬ 
ary.  However,  some  experimental  results  indicate  that  under  the  same 
load  characteristics  and  loading  sequence,  different  materials  in  some 
cases  do  not  affect  Q  significantly  [2].  It  means  that  In  some  cases 
different  materials  do  not  affect  this  analogy  method  significantly. 
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AN  ANALOGY  METHOD  FOR  CRACK 
INITIATION  LIFE  PREDICTION 

Zhang  Fuze 

(Air  Force  Research  Iexlituie) 

Abstract 

A  method  for  crack  initiation  life  prediction  of  components  is  presented. 
This  method  employs  the  life  of  a  given  component  obtained  from  endurance 
tests  under  spectrum  load  to  predict  the  life  of  another  untested  same  type 
component  by  analogy. 

The  formula  for  analogy  life  prediction  is  deduced  from  Miner's  theory 
and  formula  as  followsi 


where  */,  ■£«,  A  <7,  and  X  are  load  cyclic  number,  stress  amplitude,  overload  incre¬ 
ment  and  life  cyclic  number  of  the  given  component  in  i  level;  n S'„,  Ay*  and 
X'  are  those  of  the  unknown  component;  and  Ay.  are  any  constant  stress 
amplitude  and  overload  increment  (they  may  be  taken  as  unit  in  calculation); 
m„mj  are  constant  (generally  between  5  and  10)  and  may  be  estimated  from 
the  S~N  curve  or  by  experiments.  In  life  prediction  with  this  formula  it  is 
unnecessary  to  choose  the  constantQin  Miner's  formula,  therefore  its  accuracy 
may  be  improved. 

The  lives  of  large  components  of  two  aeroplanes  and  the  specimens  under 
five  different  loading  conditions  are  predicted  by  this  analogy  method.  The 
results  arc  quite  consistent  with  fatigue  test  results  and  the  comparison  is  shown 
in  tables  (  4  )  and  (  3  ).  This  demonstrates  that  the  formula  presented  in  this 
paper  is  more  accurate,  particularly  for  components  analogous  to  given  ones, 
and  may  give  quite  satisfactory  result?  provided  the  value  of  the  m  is  estima¬ 
ted  appropriately.  According  to  Miner's  theory  it  is  possible  to  apply  this 
analogy  method  to  components  different  from  given  ones,  but  it  is  necessary  to 
make  appropriate  correction  in  actual  use. 
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PLOW  MECHANISM  AND  EXPERIMENTAL  INVESTIGATION  OP  ROTATING 
STALL  IN  TRANSONIC  COMPRESSORS 


Beijing  Institute  of  Aeronautics  and  Astronautics, 

Lu  Yajun  and  Zhang  Sbunlin 

ABSTRACT 

In  this  paper,  a  vortex  system  flow  model  has  been  established 
for  rotating  stall  based  on  vortex  theory.  A  comparison  of  the  calcul¬ 
ated  flow  field  using  that  flow  model  with  experimental  data  was  per¬ 
formed.  This  comparison  indicated  that  its  results  agreed  with  the 
actual  flow  condition  better  than  the  small  linear  perturbation  theory 
and  other  methods. 

This  paper  gives  a  detailed  experimental  investigation  on  the  trans¬ 
sonic  rotor  undergoing  rotating  stall  at  non-designed  rotation  speed 
(8000  revolutions/minute;  13000  revolutions/minute;  15000  revolutions/ 
minute).  Some  new  flow  phenomena  worth  further  studying  were  observed. 

FOREWORD 

Due  to  the  complexity  of  the  physical  phenomenon  of  compressor 
rotating  stall  which  is  affected  by  many  factors  including  some  factors 
on  board,  despite  the  fact  that  for  many  years  people  tried  to  determine 
the  relations  between  the  aerodynamic  parameters  and  the  geometric  para¬ 
meters  of  the  compressor  and  the  characteristic  parameters  under  the 
working  conditions  of  rotating  stall,  no  satisfactory  result  has  been 
obtained.  Recently,  along  with  the  improvement  of  the  characteristics 
of  axial  compressors,  rotating  stall  and  this  type  of  nonsteady  flow 
phenomena  have  a  more  serious  influence  on  the  compressor  characteris¬ 
tics  and  destructive  effect.  It  frequently  causes  the  blade  of  the 
compressor  to  break.  It  even  resulted  in  major  incidents  involving 
loss  of  airplanes  and  human  lives.  Therefore,  the  study  of  rotating 
stall  type  of  nonsteady  flow  phenomena  is  one  of  the  important  subjects 
in  aeronautical  technology. 


This  paper  carried  out  an  analytical  study  on  the  flow  character¬ 
istics  under  the  working  condition  of  rotating  stall.  Using  the  vor¬ 
tex  theory,  a  vortex  system  flow  model  and  the  flow  field  computation 
method  were  established  under  rotating  stall  conditions.  Based  on  this 
computation  method,  a  detailed  theoretical  calculation  of  the  flow 
fields  of  a  low  speed  rotor  and  a  transonic  rotor  were  carried  out  and 
the  computed  flow  fields  agreed  with  the  actually  measured  flow  field. 
This  indicated  that  this  flow  model  and  flow  field  computation  method 
are  reasonable.  This  paper  also  contained  an  analytical  investigation 
of  oscillograms  of  the  stall  characteristic  parameters  as  a  function  of 
time  recorded  during  the  onset,  development  and  disappearance  period  of 
rotating  stall  of  a  transonic  rotor.  The  following  new  flow  phenomena 
worth  further  investigation  are:  (1)  the  "irregular  separation"  (i.e., 
serious  separation  of  individual  blades)  phenomenon  prior  to  the  onset 
of  rotating  stall  and  (2)  periodical  variation  of  circumferential  width 
of  the  stall  cell  with  time  together  along  with  the  periodical  oscilla¬ 
tion  of  the  width  of  the  stall  cell  in  the  radial  direction  of  the  blade 
in  the  duration  in  which  the  stall  flow  state  was  varied.. 
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I.  THE  ESTABLISHMENT  OF  PHYSICAL  FLOW  MODEL  UNDER  STALL  CONDITIONS 
AND  THE  COMPUTATION  EQUATION 

In  order  to  further  understand  the  flow  characteristics  and  influ¬ 
encing  factors  under  rotating  stall  conditions,  a  flow  model  capable  of 
reflecting  this  type  of  flow  phenomenon  characteristic  must  be  estab¬ 
lished.  Based  on  the  experimental  data  obtained  under  rotating  stall 
conditions  of  the  abrupt  type,  we  know  that  under  these  stall  condi¬ 
tions,  there  exists  vortex  motion  both  up  and  down  stream  of  the  blade 
array.  Therefore,  the  flow  model  for  this  condition  founded  based  on 
vortex  theory  apparently  can  better  reflect  the  characteristics  of  this 
type  of  flow  phenomenon  than  other  methods.  In  an  axial  compressor, 
the  axial  symmetry  of  the  entire  flow  field  is  destroyed  once  the 
rotating  stall  phenomenon  is  produced.  If  the  coordinate  system  is 
fixed  in  the  rotating  separation  region,  then  the  entire  flow  field  is 
divided  into  two  completely  different  diffracted  flow  regions.  In  the 


main  flow  region  outside  the  separation  region,  the  diffracted  flow  and 
moment  of  the  flow  are  basically  the  same  as  a  certain  steady  condition. 
In  the  low  velocity  flow  in  the  separation  region,  the  strong  separa¬ 
tion  of  the  flow  made  the  diffracted  flow  loop  moment  very  small. 

Because  of  the  relative  motion  between  the  blade  array  and  the  separa¬ 
tion  cell,  therefore,  each  blade  in  the  blade  array  must  continuously 
enter  and  exit  this  separation  cell  sequentially.  When  each  blade 
enters  and  exits  the  separation  cell  once,  the  flow  diffraction  moment 
on  the  blade  must  vary  once.  This  variation  is  the  source  of  the  period¬ 
ical  oscillation  force  on  the  blade.  When  the  blade  enters  the  separa¬ 
tion  cell  from  one  side  boundary  of  the  main  flow  region,  the  diffraction 
flow  moment  accompanying  the  blade  varies  from  large  to  small.  There¬ 
fore,  it  forms  a  vortex  flowing  from  the  tail  of  the  blade  toward  down¬ 
stream  of  the  blade  array  with  diffraction  flow  moment  in  opposite 
direction  to  that  of  a  blade  in  the  main  flow  region.  This  phenomenon 
is  similar  to  that  of  a  peeling  vortex  at  the  tail  of  the  wing  during 
aircraft  landing.  When  the  blade  leaves  from  the  other  side  boundary 
of  the  separation  region,  the  accompanying  diffraction  flow  moment  on 
the  blade  varies  from  small  to  large,  then  a  vortex  begins  to  flow  from 
the  tail  whose  direction  is  opposite  to  that  of  the  diffraction  flow 
moment  on  the  blade  in  the  main  flow  region.  This  phenomenon  is  similar 
to  the  formation  of  an  adherent  vortex  on  the  airplane  wing  during  '■ake- 
of'  and  simultaneously  a  take-off  vortex  is  left  behind  on  the  ground. 
Because  the  blades  in  the  blade  array  must  enter  and  exit  from  both 
sides  of  the  separation  regions  continuously  and  sequentially,  there¬ 
fore,  two  stable  vortex  streets  are  formed  on  either  side  of  the  separa¬ 
tion  region.  If  the  effect  of  viscosity  must  be  considered,  then  the 
Intensity  of  the  two  vortical  streets  downstream  from  the  blade  array 
will  gradually  decrease.  In  this  paper,  the  flow  model  uses  two  finite 
length  constant  intensity  vortical  streets  to  consider  the  effect  of 
viscosity.  The  relative  motion  between  the  low  velocity  flow  (or  even 
reverse  flow)  in  the  separation  region  and  the  incoming  flow  from  In 
front  of  the  blade  array  will  cause  the  existance  of  vortical  motion  in 
the  upstream  of  the  blade  array  (in  the  upstream  vortical  street).  The 
shape  of  this  vortical  street  is  related  to  the  size  of  the  separation 
region  and  the  flow  velocity  in  the  same  region.  In  order  to  simplify 
the  flow  model  and  the  complexity  of  the  computation,  let  us  use  a  finite 
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extension  of  the  downstream  vortical  street  towards  upstream  to  consi¬ 
der  the  effect  of  vortex  in  the  upstream  region.  Figure  1  gives  the 
vortex  system  flow  model  established  in  this  work.  From  this  vortex 
flow  model,  it  is  possible  to  deduce  the  following  flow  computation 
equation: 


when 
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in  these  equations,  T — Intensity  of  point  vortex 

t — cascade  distance  of  blade  array 
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Fig.l  A  flow  model  of  vortex  system  in  rotating  stall  condition 

Key:  1 — T  vortex;  2 — r  blade 

The  subscripts  u  and  a  represent  the  components  on  the  x  and  y-axis, 
respectively.  Similar  symbols  apply  to  whatever  follows: 
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c^ — x  coordinate  of  a  blade  in  the  separation  region 
L — circumference  of  the  blade  cascade  loop 
L30 — nominal  width  of  the  separation  region 

^•“^jCOsC^  +  S)  +Fscos  5 
fr,.-£/1sin(*-&)  +  risin  S 

2 

n  m  -  oo 
v  00 

p-  -  I  Yl  y  (  sina^aSinaC.  \ 

*  2*  ,  m-aa'  S»>a*«S“a»-  J 

6 — azimuth  angle  of  the  vortex  street 
y — strength  of  the  vortex  street 


9*-.  0»»,  <*£•,  oj,  and  <*!•  are  the  geometrical  parameters  in  the 
vortex  system  model  as  shown  in  Figure  1. 
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where  v — flow  velocity  in  the  main  flow  region 

Ca — average  value  of  gas  inlet  axial  velocity  of  the  rotor 
na — number  of  computational  points  on  the  cross-section  of 
computation 

nu — number  of  computational  points  on  the  cross-section  of  compu¬ 
tation  outside  the  separation  region 


The  detailed  computational  procedure  regarding  the  two-dimensional 
flow  field  is  reported  in  [1].  The  velocity  field  of  a  transonic  rotor 
computed  based  on  the  above  method  is  shown  in  Figure  2.  In  order  to 
verify  the  accuracy  and  adequacy  of  the  model,  its  actual  measured 
velocity  field  is  given  in  Figure  3.  Through  the  comparison  and  anal¬ 
ysis  of  the  velocity  field,  it  can  be  concluded  that  this  model  is  rea¬ 
sonable. 
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Fig. 2  The  calculated  flow  field  is  transonic  rotor  in  routing  stall  at  •  -  13000rp«n  with 

flow  coefficient  Cmm  0.263 
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Fig. 3  The  measured  flow  field  is  transonic  rotor  in  routing  stall  at  « *  13000rpm  with  ma 

flow  coefficient  Ce*  0.263 


II.  EXPERIMENTAL  APPARATUS  AND  MEASURING  EQUIPMENT 


The  test  stand  is  a  single  or  dual  level  compressor  test  equipment. 
The  test  stand  has  a  given  power  of  about  2000  horsepower.  The  given 
rotating  speed  is  22000  rpm.  The  experimental  flow  range  is  0-20  kg/sec. 
The  transonic  test  compressor  has  a  designed  rotating  speed  of 
22000  rpm  and  flow  rate  at  13.6  kg/sec.  The  extreme  pressure  ratio 
under  design  condition  is  ^1.6.  The  relative  Mach  numbers  at  the  tip 
and  the  root  of  the  blade  are  1.4  and  0.9  respectively. 

In  order  to  measure  the  steady  state  parameters,  such  as  the  total 
static  pressure  at  the  inlet  and  outlet  of  the  rotor  and  the  flow  direct¬ 
ion  at  the  outlet  of  the  rotor,  a  six-point  pressure  probe  was  installed 
on  the  cross-section  of  the  inlet  and  the  outlet.  In  addition,  four 
static  pressure  holes  were  placed  on  the  inner  and  outer  walls  at  the 
cross-section  of  the  inlet  and  outlet.  A  single-point  total  static 
pressure  combination  probe  which  was  capable  of  tracing  the  direction 
of  the  flow  was  installed  on  the  cross-section  of  the  outlet  of  the  rotor. 
The  flow  of  the  rotor  was  measured  through  the  inlet  flow  tube.  The 
steady  state  parameters  were  acquired  and  recorded  by  a  model  XJ-100  auto¬ 
matic  cycling  measuring  instrument.  Some  of  the  low  pressure  parameters 
were  obtained  directly  from  water  displacement. 

Dynamic  parameters  under  stall  conditions — the  measurements  of  flow 
velocity  and  total  gas  pressure,  were  done  by  using  model  55M  hot  fila¬ 
ment  wind  velocity  meter  of  the  DISA/CTA  series  and  the  Chinese  LDY  6-4 
pressure  pulser  as  its  sensor.  The  dynamic  parameters  were  acquired 
and  recorded  using  the  TEAC/R510  magnetic  data  recorder  for  the  entire 
duration.  During  the  whole  experiment,  a  RS-1  oscilloscope  was  used  to 
observe  and  monitor  the  dynamic  parameters  on  site.  The  measurements 
of  dynamic  parameters  and  the  treatment  of  data  were  reported  in  [2]. 

III.  EXPERIMENTAL  RESULTS  AND  ANALYSIS 
* 

1.  The  actual  measured  flow  field  and  its  analysis  of  a 
transonic  rotor  under  rotating  stall  condition 
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The  actual  measured  basic  element  level  flow  field  at  the  tip  of 
the  blade  for  a  transonic  rotor  is  shown  in  Figure  3  under  the  rotating 
stall  condition.  This  flow  field  was  measured  when  the  transonic  rotor 
was  at  13000  rpm  with  a  flow  coefficient  c.»0.263  .  From  the  flow  field 
diagram,  it  is  clearly  shown  that  there  are  reverse  flow  regions  both 
on  the  front  and  rear  cross-sections  of  the  rotor.  The  existence  of 
the  reverse  flow  areas  would  drive  the  downstream  gas  which  had  already 
been  heated  and  pressured  by  the  rotor  back  to  the  upstream  region. 


This  is  the  major  reason  for  the  melting  of  the  rotor  blade.  The  higher 
the  pressure  increase  ratio  is,  the  shorter  the  time  required  to  melt 
the  blade  due  to  the  reversed  i*low.  It  is  even  possible  to  create 
this  type  of  destructive  phenomenon  in  several  tens  of  seconds.  For 
a  rotor  with  axial  gas  inlets,  when  reverse  flow  occurs,  there  exists 
vortex  motion  in  the  upstream  region  of  the  rotor.  By  comparing  this 
flow  field  diagram  with  that  reported  in  [3]  for  an  actual  measured  flow 
field  of  a  low  velocity  rotor,  it  can  be  concluded  that  the  downstream 
reverse  flow  region  of  the  transonic  rotor  dissipates  faster  than  that 
of  a  low  velocity  rotor.  In  the  flow  field  of  a  transonic  rotor,  the 
reverse  flow  region  downstream  from  the  rotor  disappears  very  rapidly. 

On  the  N  cross-section  170  mm  from  the  rotor,  reverse  flow  cell  ceased 
to  exist.  For  low  velocity  rotors  the  downstream  reverse  flow  region 
shrinks  slowly.  It  extends  to  the  exit  channel  without  dissipating  com¬ 
pletely.  This  situation  indicates  that  the  vortex  intensity  of  the 
transonic  rotor  decreases  more  drastically  than  that  of  the  low  velo¬ 
city  rotor.  Under  rotating  stall  conditions  of  the  transonic  rotor, 
the  fact  such  that  vortical  motion  exists  both  upstream  and  downstream 
from  the  rotor  and  the  vortex  intensity  downstream  from  the  rotor 
decreases  rapidly  has  significant  meaning  to  the  correction  and  improve¬ 
ment  of  the  theoretical  flow  model. 

2.  The  process  and  analysis  of  the  onset,  development  and 
disappearance  of  rotating  stalT 

(1)  Experimental  results  and  analysis  at  8000  rpm 

The  pressure  signal  oscillograms  from  steady  state  to  the  transi¬ 
tion  period  of  stall  condition  are  shown  in  Figure  4.  These  probes 
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F.g,4  The  pressure  sigwl  osciUo,r«as  »t  »-8000rpm  in  stall  condition  and  in  blade  irregular 

separation  condition 

Key:  1  1 — time  signal;  2 — pressure  signal  (1);  3 — pressure  signal  (2); 

4 — flow  signal;  5 — rotating  speed  signal 

were  installed  at  the  blade  tip  about  10  mm  away  from  the  outer  wheel 
shell.  From  this  figure,  it  can  be  observed  that  in  the  transition 
period  from  steady  state  to  stall  condition  the  average  value  of  the 
pressure  signal  remained  unchanged.  However,  in  the  transition  period, 
the  pulse  amplitude  of  the  pressure  signal  increased  significantly. 

This  flow  phenomenon  which  showed  basically  the  same  average  pressure 
before  going  into  rotating  stall  as  in  steady  state  together  with  signi¬ 
ficant  increase  in  pulse  amplitude  was  called  "irregular  separation" 
in  [2].  The  enlarged  "irregular  separation"  pressure  oscillograms  are 
shown  in  Figure  5.  From  this  figure,  we  found  that  the  pressure  oscill¬ 
ogram  contained  17  pulses  in  the  period  of  each  revolution.  The  17 
pulses  have  the  same  cycle  period  as  that  of  the  rotor.  Because  there 
were  17  blades  on  the  rotor,  therefore,  it  is  convincing  to  say  that  66 
"irregular  separation"  is  basically  a  flow  phenomenon  of  severe  separa¬ 
tion  on  each  blade  before  the  rotor  enters  rotating  stall  conditions. 
Because  the  blade  at  the  tip  is  thin  and  curved  in  shape  for  a  transonic 
rotor,  this  further  explained  the  conclusion  obtained  based  on  low  velo¬ 
city  rotor  studies  that  a  relatively  long  period  of  "irregular  separation" 
is  frequently  found  for  basic  elements  with  small  curvature  angles. 

The  same  conclusion  applies  to  transonic  rotors. 
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Figure  5.  The  pressure 
signal  oscillograms  at 
n  *  8000  rpm  in  blade 
irregular  separation 
conditions 

Key:  1 — time  signal;  2 — pressure 
signal  (1);  3 — pressure 
signal  (2);  4 — flow  signal; 
5 — rotating  speed  signal 


Prom  the  "irregular  separation"  condition  to  the  rotating  stall 
condition,  the  average  value  of  the  pressure  signal  decreased  signifi¬ 
cantly  and  showed  a  two  zone  rotating  stall  phenomenon.  Its  rotating 

velocity  3separatIm>  Hieparation/%tor  13  °-74-  Continuing  save  the 
flow,  the  separation  region  changes  from  2  to  1J  however,  ^separation 
remained  to  be  0.74.  This  condition  continued  until  the  exhaust  flow 
saving  valve  was  completely  shut.  At  this  rotating  speed,  the  compressor 
never  had  any  coughing  vibration. 

(2)  Experimental  results  and  analysis  at  13000  rpm 

Two  pressure  probes  were  also  placed  at  the  tip  of  the  blade  in  a 
similar  manner.  When  the  flow  valve  was  closed  to  enter  the  stall  trans¬ 
ition  period,  the  longitudinal  pulse  amplitude  of  the  pressure  oscillo¬ 
grams  began  to  increase  significantly  to  produce  the  "irregular  separa¬ 
tion"  phenomenon  of  severe  separation  on  each  blade  (see  the  first  seg¬ 
ment  of  the  oscillograms  in  Figure  6).  Immediately  before  rotating 
stall  occurred,  a  new  flow  phenomenon  involving  the  periodical  variation 
of  the  average  value  of  the  pulse  amplitude  existed  (see  the  second  seg¬ 
ment  of  the  oscillograms  in  Figure  6).  Because  its  variation  period 
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Fig, 6  The  pressure  signal  oscillograms  st  nm  13000rpm  in  blade  irregular  separation  condition  and 

in  rotating  stall  condition 

1—  Pressure  signal  oscillogram  in  blade  irregular  separation  condition  with  similar  mean  ralue  of 
amplitude* 

2 —  Pressure  signal  oscillogram  in  blade  irregular  separation  with  amplitude  mean  value  earring 
periodically! 

3 —  Pressure  signal  oscillogram  in  rotating  stall  condition.  r 


Key:  1 — time  signal;  2 — pressure  signal  (1);  3 — pressure  signal  (2); 

4 — flow  signal;  5 — rotating  speed  signal 

coincided  with  that  of  the  rotation,  therefore,  it  was  not  rotating 
stall.  It  still  was  severe  separation  in  the  channels  of  each  blade. 
Furthermore,  it  must  be  explained  that  the  flow  of  gas  through  the 
channels  of  each  blade  is  no  longer  the  same.  The  appearance  of  this 
phenomenon  is  the  sign  that  rotating  stall  condition  will  soon  happen 
(see  the  transition  between  the  second  and  the  third  segments  of  the 
oscillograms  in  Figure  6).  This  flow  phenomenon  which  occurs  immed¬ 
iately  prior  to  rotating  stall  h&s  never  been  recorded  in  the  literature. 
Its  pressure  variation  oscillogram  is  shown  in  Figure  6. 

Continuing  to  cut  the  flow-  since  the  appearance  of  the  "irregular 
separation"  phenomenon,  it  then  went  into  a  steady  single  zone  stall 
condition,  an  alternating  single  and  double  zone  variable  condition  and 
a  steady  double  zone  condition  sequentially.  In  the  process  of  leaving 
the  stall  condition,  the  transition  from  double  zone  to  single  zone 
condition  and  the  steady  single  zone  condition  appeared  in  that  order. 
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Fig.  7  The  asciilogranu  measured  by  two  hot-wire  anemometer  probes  ia  transient  process  of 

rotating  stall  at  »“13000rpm 

The  npper  one  — -  the  probe  located  at  the  root  of  the  blade* 

The  lower  one— “the  probe  located  at  the  tip  of  the  blade. 

1 — time  signal;  2 — hot-wire  anemometer  signal  (1);  3 — hot-wire 
anemometer  signal  (2);  4 — flow  signal;  5 — rotating  speed  signal 
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Fig. 8  The  oecillograms  measured  by  two  pressure  transducers  in  transient  process  of  rotating  stall  at 
a*13000rpm  (transducers  located  at  the  tip  of  the  blade) 

1 — time  signal;  2 — pressure  signal  (1);  3 — pressure  signal  (2); 
4 — flow  signal;  5 — rotating  speed  signal 


Finally,  it  excited  the  steady  single  zone  stall  condition  to  a  steady 
state  condition.  In  the  stall  condition  described  above,  the  rotating 
speed  lkeparati0n  in  the  separation  region  was  about  0.7. 

In  the  transition  state  under  which  the  stall  flow  state  varied, 
the  oscillograms  measured  by  the  hot-wire  anemometer  probes  (Figure  7) 
and  the  two  pressure  transducers  (Figure  8)  were  recorded.  The  two  hot¬ 
wire  anemometer  probes  used  to  measure  the  two  oscillograms  in  Figure  7 
were  installed  at  the  blade  tip  and  blade  root  of  the  rotor  respectively 


The  angle  between  them  is  90°.  The  two  pressure  transducers  used  to 
measure  the  two  oscillograms  in  Figure  8  were  installed  at  the  roots 
of  the  blades.  Their  angle  was  also  90°.  The  four  probes  were  distri¬ 
buted  on  the  same  measurement  cross-section.  From  Figure  8,  it  is 
clear  that  the  two  oscillograms  are  basically  the  same  in  shape.  There 
was  only  a  small  phase  difference  in  the  transverse  position  of  the 
oscillogram.  It  is  apparent  that  the  transverse  distance  was  caused  by 
the  90°  angle  between  the  two  probes.  Now  we  analyze  the  variation 
of  the  two  anemometer  probe  oscillograms.  The  widths  of  the  separation 
regions  at  the  tip  and  the  root  of  the  blade  were  varied  periodically. 
Their  variation  period  is  much  higher  than  that  of  the  rotor  rotation. 

If  the  two  oscillograms  were  moved  laterally  by  half  the  period  of  the 
variation  of  the  width  of  a  separation  region,  we  can  see  that  when 
the  width  of  the  separation  region  at  the  root  changed  from  narrow  to 
wide,  the  width  of  the  separation  region  at  the  tip  changed  from  wide 
to  narrow  and  the  velocity  average  from  small  to  large.  This  is  true  for 
vice  versa.  Thus,  when  the  width  of  the  separation  region  at  the  root 
reached  the  widest  and  narrowest  values,  the  separation  region  width  at 
the  blade  tip  happened  to  be  the  narrowest  and  the  widest,  respectively. 
Combining  the  two  oscillograms  obtained  using  two  hot-wire  anemometer 
probes  at  different  blade  heights  In  Figure  7  and  the  two  oscillograms 
obtained  using  the  two  pressure  transducers  placed  at  the  tip  in  Figure 
8  for  analysis,  it  is  clear  that  in  the  transition  state  the  width  of  the 
separation  region  not  only  had  periodical  variation  along  the  circumfer¬ 
ential  direction  but  also  had  periodical  vibration  along  the  radial 
direction  of  the  blade.  The  circumferential  and  radial  frequencies 
are  the  same  at  about  18  Hz.  As  long  as  the  stall  state  begins  to 
change,  it  is  possible  to  create  this  effect.  It  appears  that  this 
phenomenon  indicates  that  under  rotating  stall  conditions,  although  the 
total  flow  across  a  certain  cross-section  of  the  rotor  is  invariant, 
yet  the  flow  of  elementsmay  change  with  time.  The  above  flow  phenomenon 
which  occurred  In  transonic  rotors  with  larger  blade  twist  along  the 
blade  height  direction  has  not  been  reported  in  the  literature  both  in 
China  as  well  as  abroad.  Its  discovery  and  further  verification  will 
provide  a  new  and  valuable  stimulation  in  the  area  of  blade  stress 
analysis.  It  also  gave  a  new  understanding  of  unsteady  flow  structure. 
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Fig. 9  The  oscillograms  measured  by  tiro  hot-wire  aaemometer  probes  ia  transient  process  of  ratntinp 

stall  at  *alSOOOrpm  - 

The  upper  one-- the  probe  located  at  the  middle  of  the  blade i 
The  lower  one  . the  probe  located  at  the  root  of  the  blade. 

Key:  1 — time  signal;  2 — hot-wire  anemometer  probe  signal  (1); 

3 — hot-wire  anemometer  probe  signal  (2);  4 — flow  signal; 

5 — rotating  speed  signal 

(3)  Experimental  results  and  analysis  at  15000  rpm 

In  the  15000  rpm  experiments,  with  the  exception  of  placing  a  hot¬ 
wire  anemometer  probe  at  the  center  of  the  blade  Instead  of  at  the  top 
of  the  blade,  the  remaining  measuring  positions  were  unchanged.  The 
entire  experimental  process  and  the  conclusions  obtained  from  the  tests 
are  similar  to  those  of  13000  rpm.  Due  to  the  space  limitation,  they 
will  not  be  described  in  detail.  In  this  section,  only  the  velocity 
oscillograms  obtained  from  two  hot-wire  anemometer  probes  in  the  stall 
state  variation  process  are  given  (see  Figure  9). 

CONCLUSIONS 

Through  the  study  of  the  transonic  rotor  rotating  stall  phenomenon 
in  this  paper,  the  following  conclusions  can  be  obtained: 


1.  The  vortex  system  flow  model  established  on  the  basis  of  the 
vortex  theory  is  reasonable  for  rotating  stall  conditions.  The  two 
dimensional  flow  fields  obtained  based  on  model  computation  agree  with 
the  measured  flow  field. 


2.  Prom  the  measured  flow  field  obtained  experimentally  under 
rotating  stall  condtiions  for  transonic  rotors,  it  is  clear  that  re^ 
verse  flow  region  exists  both  upstream  and  downstream  from  the  rotor. 

This  phenomenon  indicates  that  not  only  is  there  vortical  motion  down¬ 
stream  from  the  rotor,  but  also  upstream  from  the  rotor.  The  down¬ 
stream  reverse  flow  region  of  the  transonic  rotor  disappears  faster 
than  that  of  the  low  speed  rotor.  This  indicates  that  the  dissipation 
of  the  vortex  intensity  is  much  stronger  for  transonic  rotors  than  that 
for  low  speed  rotors. 

3.  In  the  blade  tip  area  of  a  transonic  rotor,  because  the  shape 
of  the  blade  is  flat,  straight,  pointed  and  thin,  the  "irregular  separa¬ 
tion"  phenomenon  also  occurs  in  these  basic  elements.  This  indicates 
that  the  conclusion  obtained  in  the  study  of  low  speed  rotors  that  "a 
relatively  long  period  of  "irregular  separation"  frequently  will  occur 
for  some  basic  elements  with  smaller  blade  curvature  angle  before  going 
into  rotating  stall  condition",  also  applies  to  transonic  rotors. 

k.  In  the  rotating  stall  effect  occurred  at  a  non-designed  rotat¬ 
ing  speed  for  a  transonic  rotor,  the  relative  rotating  velocity  usep_ 

aratlon  *  VparatiotAotor  r™alns  constant  which  is  about  70*  of  that 
of  the  rotor  rotating  speed.  Therefore,  the  transmitting  velocity  of 
the  blade  array  relative  to  the  separation  region  is  3055  of  that  of 
rotor  rotating  speed.  This  value  is  far  less  than  the  transmitting  velo¬ 
city  of  the  separation  region  ( 50—6056  of  rotor  rotating  speed)  for  low 
speed  rotors.  This  may  probably  be  due  to  the  effect  of  compressibility 
on  the  transmitting  velocity  of  the  blade  array  corresponding  to  the 
separation  region  which  is  significantly  increased  for  transonic  rotors. 

5.  For  a  transonic  rotor  undercertain  non-designed  rotating  speed 
and  when  the  flow  state  of  stall  varies,  the  circumferential  width  will 
vary  periodically  with  time.  Simultaneously,  it  also  oscillates  period¬ 
ically  along  the  radial  direction  with  the  flow.  Its  periodical  varia¬ 
tion  frequency  is  15-18  cycle/sec.  The  newly  observed  flow  phenomenon, 
after  further  verification,  will  provide  a  new  stimulating  source  to  the 
blade  stress  analysis  and  will  give  a  new  understanding  on  the  nonsteady 
flow  structure. 
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FLOW  MECHANISM  AND  EXPERIMENTAL 
INVESTIGATION  OF  ROTATING  STALL 

IN  TRANSONIC  COMPRESSORS 
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Lu  Yajun  and  Zhang  Shtmlin 
C Beijing  Institute  of  Aeronautics  and  Astronautics ) 

Abstract 

The  flow  characteristics  of  the  rotating  stall  is  compressors  are  studied. 
A  flow  model  for  rotating  stall  in  axial  compressors  and  its  theoretical  calcu¬ 
lation  method  are  developed  on  the  basis  of  the  vortex  theory.  By  this  method 
a  detailed  theoretical  calculation  is  completed  for  a  two-dimensional  flow  field 
in  a  transonic  rotor  at  rotating  stall.  The  flow  field  obtained  is  in  good 
agreement  with  that  measured  from  experiments.  This  shows  that  the  flow 
model  and  its  calculation  method  are  reasonable. 

The  oscillograms  of  time-varying  stall  characteristic  parameters  recorded 
in  the  onset,  growth  and  cessation  processes  of  rotating  stall  in  transonic 
ortor  have  been  analyzed  and  studied.  Some  new  flow  phenomena  deserving  of 
further  investigation  have  been  discovered  as  follows*  1)'  irregular  separation', 

j,  e.  serious  separation  of  individual  blades,  often  preceded  the  onset  of  rotating 
stall  in  compressors  with  very  small  blade-camber  angle*  2  )  periodical  vari¬ 
ation  of  the  circumferential  width  of  the  stall  cell  with  time  and  accompany¬ 
ing  periodical  oscillation  -of- the- width  of  the  stall  cell  in  the  radial  direction 
of  the  blade,  which  both  occured"ih  transient  process  of  the  rotating  stall. 
The  circumferential-  and  radiaL  oscillation  frequencies  were  the  same,  about  15 
to  18  Hz. 

The  discovered  phenomena  indicate  that  in  stress  analysis  of  the  blade  a 
corresponding  exciting  force  must  be  considered.  It  seems  that  this  discovery 
will  be  helpful  to  understanding  of  unsteady  flow  structure  in  transonic  rotors. 


A  BRIEF  DESCRIPTION  OF  THE  FIRST  AERONAUTIC  ENGINE  SYSTEM, 
STRENGTH  VIBRATION  CONFERENCE 
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The  Chinese  Aeronautics  Society  held  the  first  Engine  Structure, 
Strength  and  Vibration  Conference  in  Canton  on  December  10-15,  1981. 

There  were  98  people  from  35  units  attending.  The  conference  presented 
57  papers  with  52  papers  collected  into  a  symposium.  It  was  published 
before  the  meeting.  The  meeting  was  divided  into  structure  and  strength/ 
vibrations  for  paper  presentations  and  discussions.  This  meeting  was 
well  prepared  and  arranged.  It  paid  great  attention  to  actual  practice 
with  many  materials  and  an  academic  atmosphere.  It  should  promote 
future  research.  This  meeting  reflected  new  progress  in  the  studies 
of  structure,  strength,  and  vibration  of  engines  in  recent  years. 

The  finite  element  method  was  widely  used  in  the  study  of  engine 
structural  strength  and  vibration  analysis.  In  the  meeting,  not  only 
the  analytical  procedures  of  the  major  components  of  the  engine  were 
proposed,  but  also  some  basic  methods  applicable  to  engine  character¬ 
istics  were  presented:  for  example,  a  new  method  of  structural  auto¬ 
matic  separation,  the  finite  element  analysis  of  the  three  dimensional 
elasticity  problem,  the  turbine  elasticity  problem,  the  transition 
elements  for  analysis  of  solid-shell  combination  structure  and  tensile 
plate-plane  beam  combination  structure,  and  the  automatic  formation 
of  geometric  shapes  of  turbine  wheels,  etc. 

The  papers  in  the  areas  of  fatigue  and  crack  of  structural  compo¬ 
nents  indicated  that  research  in  this  field  is  becoming  increasingly 
important  and  a  good  foundation  was  laid.  The  papers  involved  the 
determination  of  aeronautical  engine  cyclic  load,  the  prediction  and 
experimental  evaluation  of  life  of  major  components,  and  the  analysis  of 
and  experiments  on  turbine  wheel  cracks,  etc. 

The  analysis  of  free  and  forced  vibrations  of  periodic,  on-board 
structures,  discussions  on  the  experimental  methods  to  determine 
rotor  critical  rotating  velocity,  the  analysis  of  the  combined 
wheel-blade  system  model  using  group  theory,  the  experimental  study  of 


w 


The  blade-wheel  system  coupled  vibration,  the  analysis  of  the  effect 
of  gas  inlet  irregularity  vibration  on  compressor  blade,  etc.,  repre¬ 
sented  a  new  subject  and  characteristics  of  engine  vibration  research 
in  recent  years. 

There  were  many  reports  on  the  analysis  and  experimental  study  of 
aerodynamic  characteristics  of  the  rotor-support  system  as  well. 
Especially  good  results  on  elastic  support,  compressible  hydraulic 
resistors  and  uncoordinated  motion  of  the  rotor  were  obtained. 

In  the  meeting,  many  suggestions  were  made  on  future  academic 
activities.  Many  people  believed  that  schools,  institutes  and 
factories  should  organize  to  concentrate  their  efforts  on  overcoming 
difficulties  in  order  to  produce  fruitful  results  quickly,  which 
would  also  demonstrate  the  superiority  of  the  socialist  principle 
and  save  manpower,  as  well  as  materials,  in  order  to  change  the 
poor  image  of  engines  in  the  areas  of  structure,  strength  and 
vibration.  Certain  key  problems  were  discussed  in  detail.  Finally, 
all  attendants  recommended  that  the  next  Engine  Structure,  Strength 
and  Vibration  Conference  be  held  in  1983. 
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A  KIND  OP  TRANSITION  FINITE  ELEMENT  FOR  ANALYSIS  OP  SOLID-SHELL 
OR  TENSILE  PLATE— PLANE  BEAM  COMBINED  STRUCTURE 

Nanhua  Power  Plant  Institute,  Yin  Zeyong,  Yin  Jing  and  Ren  Peizbeng 
ABSTRACT 

A  kind  of  transition  element  was  established  in  this  paper.  When 
coordinated  with  standard  isoparametric  solid  (or  plane)  elements  and 
superparametric  shell  (or  plane  beam)  elements  of  the  same  order,  it 
can  be  used  to  treat  the  non-standard  joint  combined  structure  such  as 
solid-shell  or  tensile  plate-plane  beam  structures.  The  procedure  to 
establish  this  kind  of  transition  element  was  described.  The  formula 
of  quadratic  transition  element  was  presented.  The  numerical  results 
obtained  showed  that  this  kind  of  transition  finite  elent  is  very 
effective . 

I.  INTRODUCTION 

The  engines  used  In  aeronautical  applications  have  many  Important 
parts  and  components  which  are  composed  of  two  parts  such  as  solid  and 
shell  (tensile  plate)  or  axisymmetric  solid  and  axisymmetric  shell  (axi- 
symmetric  tensile  plate).  Some  of  these  can  be  treated  as  a  tensile 
plate-plane  beam  combined  structure  on  the  same  plane  after  simplifica¬ 
tion.  Examples  are  eccentric  turbine,  convex  shoulder  blade,  etc. 

When  the  finite  element  method  is  used  for  the  analysis  of  stress  and 
vibration  of  these  structures,  if  solid  or  plane  elements  are  used  for 
the  entire  combined  structure,  then  not  only  the  data  preparation  and 
computational  load  are  increased  but  also  possibly  some  symptoms  are 
brought  into  the  equations  to  be  solved  [1]. 

If  solid  element  or  plane  elements  are  used  for  the  solid  portion 
or  tensile  plate  portion  of  combined  structures,  and  shell  elements  or 
plane  beam  elements  are  used  for  the  shell  portions  as  plane  beam  por¬ 
tion,  then  a  coordinative  distortion  problem  will  appear  between  the  two 
kinds  of  elements.  In  order  to  solve  this  coordinative  problem,  primary 
and  accessory  variable  [2,3,4]  combination  element  [5,6]  penalty 


element  [4]  and  transition  element  [7,8]  methods  can  be  used. 

The  primary  and  accessory  variable  method  and  the  penalty  method 
are  effective  means  to  process  non-standard  joint  combined  structures. 
However,  because  of  the  need  to  distinguish  the  displacement  vector  of 
the  primary  and  accessory  nodes  and  to  have  the  accessory  node  dis¬ 
placement  vector  not  to  appear  in  the  equations,  the  program  design  of 
the  primary/accessory  variable  method  is  relatively  complicated.  As 
for  the  penalty  element  method,  the  choice  of  the  numerical  value  of  the 
penalty  matrix  has  a  certain  degree  of  difficulty.  In  this  paper,  a 
kind  of  transition  element  is  established  which  is  different  from  the 
ones  described  in  [7,8].  This  kind  of  transition  element  of  any  order 
can  be  used  coordinatively  with  the  standard  isoparametric  solid  (or 
plane)  element  and  superparametric  shell  (or  plane  beam)  element  of  the 
same  order.  With  respect  to  the  program  of  the  widely  used  standard  iso 
parametric  solid  (or  plane)  element  and  superparametric  shell  (or  plane 
beam)  element,  it  is  possible  to  add  the  kind  of  transition  element 
established  in  this  paper  in  order  to  facilitate  the  analysis  of  non¬ 
standard  joint  combined  structures  such  as  solid-shell  and  tensile  plate 
plane  beam. 
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II.  ARRANGEMENT  OF  NODES,  GEOMETRIC  SHAPE  AND  DISTRIBUTION 

OF  DISPLACEMENT 

The  basic  procedure  to  establish  this  kind  of  transition  element 
can  be  described  in  the  following:  On  the  interface  (or  boundary  line) 
of  a  standard  isoparametric  solid  (or  plane)  element  (for  a  three- 
dimensional  solid  element,  sometimes  on  an  edge),  the  original  nodes  of 
the  solid  (or  plane)  element  are  eliminated  and  replaced  by  the  corres¬ 
ponding  nodes  of  a  superparametric  (or  plane  beam)  element  of  the  same 
order.  As  for  the  remaining  nodes  of  the  solid  (or  plane)  element,  the 
shape  functions  are  identical  to  the  corresponding  shape  functions  of 
standard  i&uparametric  solid  (or  plane)  element.  As  for  the  shell  (or 
plane  beam)  element  nodes,  the  shape  functions  are  the  same  as  the  cor¬ 
responding  ones  of  the  standard  parametric  shell  (or  plane  beam) 
element. 


Obviously,  the  shape  function  of  this  kind  of  transition  element 
satisfies  the  requirement  that  its  value  is  1  at  one  node  and  0  at  any 
other  node.  Its  distribution  on  the  interface  (or  boundary  line)  nodes 
of  the  superparametric  shell  (or  plane  beam)  element  distortedly  coor¬ 
dinates  with  the  standard  superparametric  element  of  the  same  order. 

As  for  the  corresponding  interface  (or  boundary  line)  which  contains 
nodes  of  an  isoparametric  solid  (or  plane)  element  alone,  it  coordinates 
distortedly  with  the  standard  isoparametric  solid  (or  plane)  element  of 
the  same  order.  Similarly,  the  two  neighborning  transition  elements 
are  also  coordinatively  distorted.  Prom  the  following  equation,  this 
coordination  can  be  verified.  In  addition,  for  this  kind  of  transition 
element,  besides  arranging  the  interface  (or  boundary  line)  of  the  shell 
(or  plane  beam)  element  nodes,  the  displacement  inside  the  element  does 
not  agree  with  the  distortion  assumption  of  the  shell  (or  plane  beam). 

In  practice,  for  solid-shell  or  tensile  plate-plane  beam  combined  struc¬ 
tures,  with  respect  to  the  shell  (or  plane  beam)  portion  near  the  Joint 
the  ideal  shell  (or  plane  beam)  distortion  assumption  should  also  not  be 
used.  Thus,  the  established  plane  transition  element  should  be  In  a 
plane  stress  or  plane  strain  state,  the  axisymmetric  transition  element 
should  be  in  axysymmetric  solid  stress  state  and  the  three-dimensional 
solid  transition  element  should  be  In  a  three-dimensional  stress  state. 
From  the  actual  equation.  It  is  possible  to  verify  that  the  displace¬ 
ment  distribution  of  this  kind  of  transition  element  can  realize  a  con¬ 
stant  strain  state. 

The  following  is  an  example  using  a  quadratic  transition  element 
to  obtain  the  expressions  of  the  geometric  shape  and  displacement  dis¬ 
tribution. 

1.  Quadratic  axisymmetric  transition  ring  element  and 
quadratic  plane  transition  element 

Here  only  the  former  Is  discussed.  The  relevant  expressions  for 
the  latter  are  identical.  The  matching  quadratic  axisymmetric  solid 
isoparametric  elements  of  the  quadratic  axisymmetric  transition  ring 
element  are  given  in  [9].  The  matching  quadratic  axisymmetric  super¬ 
parametric  elements  can  be  established  by  initiating  the  method  in  [10]. 
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However,  we  used  a  slightly  different  method  direetly  [11].  Therefore, 
the  total  coordinates  (r,z)  and  displacement  (u,v)  of  the  local  coor¬ 
dinate  point  (S,n)  in  the  axisymmetric  transition  ring  element  can  be 
determined  based  on  the  following  equations  (see  Figure  1). 


(l) 


where  r^,  z^,  ri,  z^,  u^,  v^,  u^,  v^  are  the  coordinate  components  and 
displacement  components  of  the  corresponding  nodes,  is  the  other  dis¬ 
placement  component  of  node  1.  The  definitions  of  a^,  <J>^  and  t^  are 
all  shown  in  Figure  1.  The  actual  expressions  of  the  shape  functions 
are 
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2 .  Quadratic  three-dimensional  transition  element 

The  quadratic  three-dimensional  solid  element  to  be  matched  with 
this  kind  of  transition  element  can  be  obtained  from  [12].  The  match¬ 
ing  quadratic  superparametric  shell  elements  can  be  found  in  [10]. 
Therefore,  the  total  coordinates  (x,y,z)  and  displacement  (u,v,w)  of  a 
local  coordinate  point  (£,  n,  ?)  in  this  type  of  transition  element 
can  be  determined  from  the  following  equation  (see  Figure  2): 
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III.  STIFFNESS  MATRIX,  MASS  MATRIX  AND  LOAD  VECTOR 

The  general  expression  of  the  stiffness  matrix  [k]  for  transition 
element  is  still  [12] 

C*)-  LcS)rC£»CBWr 

J  v  (6) 

where  [D]  is  the  elasticity  matrix  of  the  material  which  can  be  chosen 
based  on  the  stress  condition  of  the  transition  element.  [B]  is  the 
strain  node  displacement  matrix.  The  following  is  the  actual  express¬ 
ion  corresponding  to  the  quadratic  transition  element. 


1.  Quadratic  axisymmetrlc  transition  ring  element  and  quadratic 
plane  transition  element 


For  the  quadratic  axisymmetrlc  transition  ring  element,  here  we 
use  the  following  equations  to  define  its  strain  vector  {e}  and  node 
^placement  vector  {5}: 
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Under  the  conditions,  the  following  system  of  equations  exist: 
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where  [J]  is  the  "Jakobi"  matrix,  |J|  and  [J]-^  are  the  determinant 
and  inverse  of  [J]  respectively. 
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where  a.  y*«  *'  and  «#.  v„  u>,  are  the  coordinate-  components  and  displace- 
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vector  at  node  j ;  fj-lv./l  is  the  thickness  at  node  j  .  v^  and  v2j  are 

determined  based  on  the  following  two  equations,  respectively: 
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If  v3J  and  I  =  (1,0,0)  are  parallel,  then  j  =  (0,1,0)  is  used  to 


replace  I.  u^,  v^,  w^,  a^,  B^  are  5  displacement  components  of  node  j 
where  oij  and  Bj  are  the  rotating  angles  of  v^j  around  v?j  and  res¬ 
pectively.  The  shape  function  expressions  are: 
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If  one  wants  to  obtain  the  relevant  equations  of  the  quadratic  plane 
transition  element,  the  fourth  row  of  the  matrices  on  the  right  hand 
of  equations  (7),  (12)  and  (13)  should  be  eliminated. 

2 .  Quadratic  three-dimensional  transition  element 


Here  the  following  equations  are  used  to  define  the  strain  vector 
{e}  and  node  displacement  {6}  of  the  element: 

<e}-Ce„  e„  e„  Y^)T  (14) 
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Under  the  above  conditions,  the  following  equation  system  exists: 
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Similarly,  based  on  the  general  expression  [1],  the  mass  matrix  _ 

and  load  vector  of  the  transition  element  can  be  obtained.  It  is  not 
discussed  in  detail. 

IV .  EXAMPLES 

Due  to  the  limitation  of  space,  here  we  only  give  two  exam¬ 
ples  to  show  the  accuracy  and  application  of  this  kind  of  transition 
element.  The  analyzed  combined  structures  have  unit  thickness  and  were 
under  plane  stress  situation. 

1.  Transverse  force  exerted  on  the  end  of  staircase- 
shaped  cantilever  beam 

The  data  of  this  example  are  shown  in  Figure  3(a).  The  element 
subdivisions  are  shown  in  Figure  4(a)  and  (b)  where  (a)  used  quadratic 
plane  element,  plane  transition  element  and  plane  beam  element,  and  (b) 
used  quadratic  plane  beam  element  only.  The  deflections  of  line  AB 
calculated  using  both  subdivisions  and  its  analytical  solution  [13]  are 
all  shown  in  Figure  5. 

2 .  Plane  fork  structures  with  two  concentrated  forces 

The  data  of  this  example  are  shown  in  Figure  3(b).  Using  symmetry, 
only  half  of  the  structure  is  examined.  In  this  case  we  can  no  longer 
use  plane  beam  elements  alone  to  carry  out  the  analysis.  However,  the 
method  involving  the  use  of  transition  elements  as  shown  in  Figure  4(a) 
is  still  feasible.  The  calculated  deflections  of  line  AB  are  also 
shown  in  Figure  5.  Figure  6  shows  the  horizontal  displacements  of 
some  cross-sections. 
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A  KIND  OF  TRANSITION  FINITE  ELEMENT  FOR 
ANALYSIS  OF  SOLID-SHELL  OR  TENSILE  PLATE- 
PLANE  BEAM  COMBINED  STRUCTURES 

Yin  Zeycmg,  Yin  Jing,  Ren  Peizheng 
( Nanhua  Power  plant  Institute ) 

Abstract 

A  kind  of  transition  finite  element  is  developed.  In  coordinative  combina¬ 
tion  with  standard  iso-parametric  solid  (or  plane)  elements  and  super-parame¬ 
tric  shell  (or  plane  beam)  elements  of  the  same  order*  this  kind  of  transition 
element  can  be  used  to  deal  with  non-standard  joint  combined  structures  such 
as  solid-shell  or  tensile  plate-plane  beam  structures. 

The  procedure  to  develop  this  kind  of  transition  element  is  as  follows. 
On  an  interface  or  an  edge  of  a  standard  iso-parametric  solid  element  or  on 
a  boundary  line  of  a  standard  iso-parametric  plane  element*  the  corresponding 
nodes  of  a  super-parametric  shell  (or  plane  beam)  element  of  the  same  order 
are  placed  instead  of  the  nodes  of  the  original  solid  (or  plane)  element.  The 
shape  functions  for  remained  solid  (or  plane)  element  type  nodes  are  the  same 
as  the  corresponding  ones  for  a  standard  iso-parametric  solid  (or  plane) 
element.  As  regards  shell  (or  plane  beam)  element  type  nodes  the  same  shape 
functions  are  adopted  as  for  a  standard  super-parametric  shell  (or  plane  beam) 
element. 

The  formulation  of  equations  Tor  the  quadratic  transition  element  is  pre¬ 
sented  in  detail. 

The  numerical  results  obtained  have  shown  that  this  kind  of  transition 
finite  element  is  very  effective. 


PATTERNS  OP  MONOPULSE  ARRAYS  WITH  TRIANGULAR 
AMPLITUDE  DISTRIBUTION 

Beijing  Institute  of  Aeronautics  and  Astronautics,  Liu  Sbanwei 
A  BSTRACT 

This  paper  discuses  the  sum  and  difference  pattern  function  of  83 
monopulse  arrays  of  equally  spaced  discrete  elements.  The  excitation 
amplitudes  include  uniform,  triangular  and  V  shaped  distributions.  The 
compact  form  of  array  factor  expression  ( 4),  SC-*)  —  F(*)  — (?(*)  has  been 
obtained  by  other  people  using  either  exact  or  truncated  Z-transform. 
With  the  aid  of  equation  (8),  and  sometimes  taking  the  shifting  theorem 
into  account,  S(z)  can  be  found  simply  by  the  unilateral  Z-transform 
without  considering  whether  the  envelope  function  is  symmetrical  or 
asymmetrical.  A  numerical  example  was-  also  given.  A  study  was  carried 
out  with  respect  to  the  sum  pattern.  In  the  major  region,  the  calcul¬ 
ated  results  agreed  very  well  with  the  experimental  data.  The  differ¬ 
ence  pattern  has  also  been  verified  experimentally. 

I.  Z-TRANSFORM  AND  ARRAY  FUNCTION 

In  the  early  1960's,  D.  K.  Cheng  and  M.  T.  Ma  for  the  first  time 
successfully  applied  the  Z-transform  for  treating  the  time-space  prob¬ 
lems  to  the  pulse-data  system  in  the  analysis  of  discrete  array  patterns 
and  expanded  the  applications  of  Z-transform  theory.  Several  relevant 
papers  were  published  later  [2,3,4].  In  this  section,  the  derivation 
of  array  function  is  further  simplified. 

As  is  commonly  known,  the  linear  array  factor  of  equally  spaced 
linearly  varying  excitation  phase  radiation  elements  can  be  written  as: 

5*  ^  /<expCj*’W(oo3  0  — cos0#)3“  2  ^iZ  (1) 

»  -  •  •  ■  0 

*  “exp (—;»)  (2) 

a  -W(coe0  -co80.)  (3) 

where  (X  is  the  wavelength  in  free  space  );  d  is  the  distance 

between  neighboring  radiation  elements;  /  -*/-  l ,  0  is  the  angle  from 
the  linear  array  axis,  8  is  the  maximum  radiation  direction; 
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* 

w.- 


/  (id),  ( i -0,  1,  •••»  "  ~  1  ^  is  the  excitation  amplitude;  n  is  the 
total  number  of  radiation  elements;  /(C)  is  the  envelope  function 
of  excitation  amplitude. 

In  equation  (1)  there  are  n  terms  which  is  also  called  an  array 
polynominal.  It  can  be  rewritten  as: 

•  -  i 


where 
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$(*)-  a  * (,d)z/ 

(4a) 

i  *  0 

-F(*)-G(*) 

(4b) 

F(*)-S/(W-Z(/(t)) 

$  *0 

(5) 

*o 

a  /(*</)*”«  ZC/(C)l/(C-m/) ) 

i  «a 

(6) 

84 

l  0  ,  C  <md 

(7) 

F(z)  and  G(z)  are  the  Z-transforms  of  two  functions  and  LKC-wOis  the 
shifted  unit-step  function  of  displacement.  Because  the  exponent 
(n-1)  of  I  in  equation  (4a)  is  a  finite  number,  therefore,  S(z)  is 
called  the  finite  Z  transform  of  f (c) 

The  ordinary  Z  transform  of  f(c)  in  equation  (4b)  is  F(z)  which 
is  easy  to  locate  in  the  theory  and  tables  introduced  in  ordinary 
reference  books  [5,6].  Therefore,  the  key  to  finding  S(z)  is  to  obtain 
G(z).  It  can  be  proved  that  when  1>Q  G(z)  has  the  following  relation: 

G(*)-ZC/(t)t/(t'-nrf):)-z"ZC/(t  +*OD  (8) 

Thus,  it  is  possible  to  obtain  the  expression  G(z)  from  equation  (8). 

It  is  equal  to  the  Z-transform  of  the  function  /  (C  +«0  multiplied  by 
z”n.  This  physical  meaning  of  equation  (8)  is  obvious  which  is  the  Z- 
transform  of  the  function/ (C +»rf)af ter  a  displacement  of  nd.  Therefore, 
it  is  possible  to  conveniently  obtain  G(z)  using  the  unilateral  Z-trans¬ 
form  and  shifting  theorem. 
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II.  SUM  AND  DIFFERENCE  PATTERNS  OF  TRIANGULAR 
AMPLITUDE  EXCITATION  ARRAY 


The  distribution  of  the  triangular  amplitude  excitation  is  shown 
in  Figure  1.  Figure  1(a)  is  sum  excitation,  and  Figure  1(b)  is  differ 
ence  excitation.  The  units  are  distributed  evenly  on  the  c  axis  at  a 
distance  d.  We  assumed  the  total  number  of  element  n  is  even,  the 
array  function  is  obtained  using  equations  (4)-(8). 


,/fl  h\- 


,]r 


—  (n— l)d 


Fig .  1  Triagnler  amplitude  distribution  for  even  ■ 
C  n  )  saau  (  b  )  difference . 


The  truncated  Z  transform  of  the  left  half  of  the  array  function 
in  Figure  1  is  SL(z)  and  the  truncated  Z-transform  of  the  right  half 
from  C+(»-i)d  is  S„(z),  i.e.. 
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>-*  2  *l*+— 
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<1 

«  x  ~~2~*  1  " 

>K ']/ ( 

2  V  2 
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■*"  *  z[-t  +jltL</J“ 

Z  C-C) 
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-dx[H~2  " 

-  *  2 

djfl  2  *  2 

*  i  « 

$•.(*)+$«(*)  •$«(*)  is  the  sum  excitation  array  function; 

St(*)-S«(*)-S4(*)  is  the  difference  array  function.  By  substituting 

equations  (9)  and  (10)  into  them  and  after  rearranging,  we  get 

St(  * ) l  -*'•'*) C  l  -,*<►*>/»}/(  *  -  1  )•  ( 11 ) 

54(*  )-<<*(  l+(«  —  l)(l  —  *)**"'*— **'“*MV(  *  -  1  )*  (12) 


n-l  (xd> 


-1  ZxJ) 


1+ot-  (it- Dad 


2  Triangular  amplitude  dUtnbutioo  luperimpoaed  on  uniform  )M>t 
<«)  sum.  <b>  difference. 


Using  a  similar  method,  we  can  obtain  the  array  function  of  an 
even  number  radiation  element  linear  array  spaced  by  d  as  shown  in 
Figure  2  when  the  excitation  amplitudes  are  uniform  and  triangular 
addition. 


In  the  case  of  sum  excitation  (Figure  2a) 

y  «)-  -  [  1  ziL  j  i  -«•"'»)(  i  1 

L  *  -  1  ( 2  -  i  )*  J 


(13) 


In  the  case  of  difference  excitation  (Figure  2b) 
C  (  ,  .  ,  r  (  1  -zm/ty  ad  C(  H  -  1  )(?',/|,|-^‘*/>)t:'”1-  1  3 

5  l  2-1 


(2-  l)s 


]  (14) 


It  is  obvious  that  the  maximum  value  of  equation  (13)  occurs  at 

W')-S(') -Inf.  (  *  ) ' -§Hr ' -  *  [  1  +  (15) 


where  and  Q2  represent  the  numerator  and  denominator  of  S^.  (z)  res 
pectively.  Q1”  and  Q2"  are  the  second  order  derivatives  of  Q1  and  Q,,. 
When  z  *  1,  equation  (14)  is  equal  to  zero,  i.e.,  S4(l)—  0.  . 


For  ease  of  application,  it  is  customary  to  further  simplify  S(z) 
and  the  2  *exp (-;'«)  in  equation  (2)  is  substituted  and  restored  as  a 
function  of  u.  Through  the  simplification  of  equation  (13),  and  omitt 
ing  the  common  factor  term  which  represents  the  phase  center 

position  of  the  far  away  region  of  the  array,  the  expression  of  sum 
pattern  S(u)  can  be  obtained: 


S(«)’ 
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ad  sin - - - u 

4 _ 
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sin— 


(16) 


(17) 


If  the  array  is  formed  by  the  same  radiation  elements,  the  actual 
array  pattern  is  the  product  of  the  above  equation  and  the  element  pat 
tern.  For  example,  the  array  corresponding  to  equation  (16)  is  formed 
by  half  wave  oscillator  coaxial  array  with  a  pattern: 
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Let  us  choose  n  ■  14,  ad  ■  g-,  then  the  array  sum  pattern  becomes: 
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The  maximum  radiation  value  Is 

S.«.-S«„-S(0)-21 


(20) 


In  Figure  3  we  obtained  the  calculated  results  from  equation  (19) 
and  the  experimental  results.  The  dotted  line  In  the  figure  Is  based 
on  calculated  results.  The  solid  line  is  the  pattern  obtained  from 
the  x  band  rectangular  wave  width  longitudinal  half  wave  crack  array. 
Its  excitation  form  is  as  shown  in  Figure  2(a).  The  comparison  of  cal¬ 
culated  and  measured  results  showed  that  in  the  major  region  the  two 
basically  coincide.  The  half  power  wave  band  width  is  very  close 
which  is  about  2A0«*5.3i*  •  In  the  visible  region  of  0  -  0  —  «  (I  ±  »  K251.82' ) 
the  calculated  zero  point,  side  lobe  position  and  relative  voltage  are 
as  follows: 


|3  point  Mi,*.  ±30.92,  ±51.43,  ±78.94,  ±  102.86,  ±  129.28,  ±154.29. 
position 

Wide  lobe  ±180.0(7,  ±205.71,  ±230.73,  ±251.82, 

position  u;«»±40,  ±64,  ±90,  ±115,  ±141,  ±167,  ±192,  ±217,  ±241, 

relation  volt- ^,  Q  097 5 x  +o. 08944,  -0.05087,  +0.04599,  -0.03437, 
age  of  the  side  lobes 

+0.02992,  -0,02485,  +0.01871,  -0.01172, 


According  to  equation  (14),  after  simplification  we  can  obtain  the 
difference  pattern  of  the  excitation  array  with  a  shape  as  in  Figure 
2(b)  using  u  as  the  independent  variable: 
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Still  using  n  ■  14  and  ad  *  |  to  substitute  into  equation  (24),  we  get 
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D,(0)-limI>,(«O=  0 
u  ♦  0 


(23) 


The  results  calculated  from  equation  (22)  are  plotted  in  Figure  4 
Only  the  region  where  is  drawn.  The  region  in  which  «<0  appears 

to  be  antisymmetric.  The  validity  of  the  equations  have  been  verified. 


Fig. 3  Comparison  of  m.asured  result  with  calculated  cunre  on  radiation  pattern  at 
9375MHz  of  id-element  linear  array 


III.  SUM  PATTERNS  OF  INVERTED  TRIANGLE/OR  V  SHAPED 
AMPLITUDE  EXCITATION  ARRAY 

Based  on  the  method  Introduced  in  the  first  section,  we  can  derive 
the  array  function  of  the  inverted  triangular  excitation  as  shown  in 
Figure  5,  i.e.,  the  sum  pattern. 


When  n  is  an  odd  number 

S.(z)-adz[2-^-(z  +*-,)  +  2*‘-u'*]/(*-l)t  (24) 


When  n  is  an  even  number 

S.(  *  *  +*-)-JLyX(  1  +*—)  +  <  1  +  *  •'*)/(  z  -  1  )*  ( 25 ) 

Changing  the  variable  to  u,  and  omitting  the  terms  with  the  absolute 
value  of  the  common  exponent  equal  to  1,  equations  (24)  and  (25)  then 
are  transformed  into: 
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It  is  obvious  that  when  u  *  0,  the  radiation  has  its  maximum 

S.^=S.(0)-limS.(a)— 2SL(«“—  1  )  (28) 

tl  0  4 

S S •(.  0  )  * a  )  *“7  **  (29) 

if  •*>  o 

when  .  Naturally,  the  same  results  can  be 

obtained  using  equations  (24)  and  (25). 

The  difference  excitation  corresponding  to  the  V-shaped  sum  excita 

tion  as  shown  in  Figure  5  has  a  linear  amplitude  distribution.  The 

array  function  (i.e.  difference  pattern)  obtained  agrees  with  that  in 

[4]  and  will  not  be  repeated  here. 
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PATTERNS  OF  MONOPULSE  ARRAYS  WITH 
TRIANGULAR  AMPLITUDE  DISTRIBUTION 

Lii  Shamuei 

( Beijing  Institute  of  Aeronautics  and  Astronautics) 

Abstract 

This  paper  describes  the  sum  and  difference  patterns  of  monopulse  linear 
arrays  of  equally  spaced  discrete  inphase  elements  which  are  excited  in  such 
a  way  that  the  envelope  of  the  excitation  amplitudes  is  uniform,  triangular, 
or  V  in  shape.  Analysis  of  the  radiation  patterns  is  carried  out  by  using  exact 
or  truncated  Z-transform,  and  then  the  compact  formas  S  (e)~  (?(*), 

’  *■  *•  (  O*  is  obtained.  This  method  was  first  developed  by  D.  K.  Cheng  and 

M.  T.  Ma  in  1960.  With  the  aid  of  Eq.  (  8  )  or  the  formula  of  G  (?)  derived 
from  the  shifting  theorem*  S  (a)  can  be  found  simply  by  the  unilateral  Z- 
transform  without  considering  whether  the  envelope  function  is  symmetrical  or 
unsymmetricaL  Also  a  numerical  example  i*  given*  showing  that  the  theorem 
tically  calculated  sum  pattern  compares  quite  favorably  with  that  obtained 
from  experiments*  especially  in  the  major  region  the  test  data  and  the  ana— 
lytieal  results  have  been  found  in  fairly  satisfactory  agreement  The  difference 
pattern  also  has  been  checked  by  experiments.  The  method  proposed  here  is 
simple  and  easy  to  carry  out 
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MODIFICATION  OF  THE  WILLENBURG  AND  MAARSE  MODELS  AND 

APPLICATION  TO  THE  LIFE  PREDICTION  FOR  T1-6A1-4V  ALLOY  89 


Beijing  Aeronautical  Materials  Research  Institute, 

Zhang  Yongkui,  Gu  Mingda  and  Yan  Minggo 

ABSTRACT 

The  accuracy  of  the  overload  retardation  effect  on  fatigue  crack 
growth  is  closely  related  to  the  formula  for  the  plastic  zone  chosen 
in  the  retardation  model.  Hence,  this  paper  proposed  a  corrected  plas¬ 
tic  zone  formula  to  replace  those  in  the  Willenborg  and  Maarse  models. 

In  addition,  the  material  constants  C#  and  n*  processed  by  the  AKef.f 
treatment  using  the  Maarse  model  were  replaced  by  the  constant  ampli¬ 
tude  material  constants  C  and  n  in  the  computation.  Through  the  calcu¬ 
lation  of  the  overload  retardation  effect  of  T1-6A1-4V  titanium  plate, 
the  improved  model  shown  in  this  paper  demonstrated  that  it  not 
only  has  significantly  higher  accuracy  than  the  original  models  but  also 
is  convenient  to  use. 

I .  INTRODUCTION 

The  calculation  of  overload  retardation  effect  in  fatigue  crack 
growth  is  one  of  the  important  problems  of  crack  growth  life  prediction 
under  spectrum  loads.  Hence,  we  have  carried  out  experimental  studies 
on  overload  retardation  characteristics  of  T1-6A1-4V  titanium  plate 
under  various  overload  ratios  and  various  crack  lengths.  In  addition, 
the  Wheeler  model,  Willenborg  model,  Matsuoka  model  and  Maarse  model 
were  chosen  to  calculate  retardation  effect  quantitatively.  Comparisons 
and  evaluation  of  the  above  four  retardation  models  were  also  carried 
out  [1].  The  study  results  indicate: 

(1)  the  fatigue  crack  growth  form  of  the  titanium  plate  under  test 
was  primarily  plane  strain  or  a  mixed  model  primarily  with  plane  strain; 

(2)  the  accuracy  of  the  Maarse  model  and  Willenborg  was  more 
inferior; 

(3)  the  selection  of  the  formula  in  the  plastic  zone  has  a  very 
significant  effect  on  the  estimation  of  the  retardation  effect. 
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As  is  well  known,  the  overload  plastic  zone  is  one  of  the  impor 
tant  factors  controlling  the  retardation  effect.  The  commonly  used 
formula  for  the  dimension  of  the  plastic  zone  is  [2]: 


where  a  is  a  constant:  Under  plane  stress  condition,  ot  —  2 ;  under 
plane  strain  condition,  a  =  6 . 

In  reality,  the  specimens  seldom  were  found  in  ideal  plane 
stress  or  plane  strain  conditions.  More  of  them  are  situated  in 
between — mixed  type  stress  state.  Especially  for  specimens  or  struc¬ 
tures  during  service,  along  with  the  growth  of  the  crack,  the  stress 
condition  varied  continuously  from  plane  strain  -*•  mixed  type  -*•  plane 
stress.  Therefore,  the  use  of  a  fixed  constant  a  to  calculate  the 
dimension  of  the  plastic  zone  R  under  varying  stress  conditions  would 
definitely  cause  larger  errors. 

Received  June  1981. 
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Hence,  in  order  to  reflect  the  characteristics  of  the  fatigue 
growth  process  and  the  continuity  of  the  stress  state  transformation 
process,  and  to  further  improve  the  prediction  accuracy  of  the  retard¬ 
ation  model,  we  made  attempts  to  improve  the  Willenborg  and  Maarse  models 
by  correcting  the  plastic  zone  dimension  formula. 

II.  IMPROVEMENT  FOR  WILLENBORG  RETARDATION  MODEL 

In  the  Willenborg  mciel  [3]»  the  effective  stress  intensity  factor 
region  AKeff.  and  effective  stress  ratio  AReff,  were  adopted.  Because  of 
ease  of  use,  it  has  received  wide  attention.  However,  tc  further 
improve  the  accuracy  of  the  Willenborg  model  for  the  prediction  of  over¬ 
load  retardation  effect,  several  researchers  used  different  angles  to 
carry  out  the  improvement  work  and  obtained  relatively  more  satisfactory 
results  [4-6].  The  improvement  of  the  Willenborg  model  in  this  paper  is 
limited  to  the  corrective  computation  of  the  constant  a  in  the  adopted 
plastic  zone  dimension  formula  (1). 
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As  is  well  known,  for  a  given  material  the  change  of  stress  condi¬ 
tion  during  crack  growth  is  mainly  dependent  upon  the  variations  of  the 
stress  intensity  factor  region  AK  and  the  stress  ratio  R.  Hence,  when 
AK  is  very  small  or  approaching  the  threshold  value  AKfch,  it  can  be  con¬ 
sidered  in  a  plane  strain  state  and  the  a  in  equation  (1)  should  be  6. 
When  AK  is  very  large  or  very  close  to  AKQ  (i.e.,  Kmax  is  very  close  to 
Kq),  it  can  be  considered  in  a  plane  stress  state  and  the  a  should  be  2' 
in  equation  (1).  The  intermediate  section  is  in  a  mixed  stress  state 
and  the  a  should  be  a  fixed  function  of  AK.  Therefore,  we  can  assume 
that  when  the  stress  relaxation  condition  is  not  considered 


1  +2  S 


where  K  is  the  fracture  modulus, 
o 


~  AK-AK,k  _  ( 1  -R)Km^-AK,k 

A  K.  (l+i?)A. 


A' _ —  - 


AK, 


( 1  —  R ) 


A'. 


(2) 


(3) 


where  Kq  is  the  fracture  modulus . 

Equations  (1),  (2)  and  (3)  are  the  improved  plastic  zone  calcula¬ 
tion  formula  of  Willenborg  model.  Obviously,  a  is  no  longer  a  fixed 
constant;  it  is  a  function  of  Kmax  and  R. 

The  physical  meaning  of  S  can  be  envisioned  as  the  ratio  of  plane 
stress.  Prom  equations  (2)  and  (3),  we  get: 

under  plane  strain,  AK  is  small  or  AK  -*■  AKth, 

S  -*■  o  ,  a  +  6 

under  plane  stress,  AK  is  large  or  Kmax  -*■  Kq, 

S  -*■  1 ,  a  •*  2 

In  mixed  stress  conditions,  o<S<l,  2<a<6. 


In  order  to  verify  the  validity  of  equation  (3),  an  example  is 
given  here.  The  experimental  data  were  taken  from  [7];  the  material 
used  was  3.2  mm  thick  2024-T3  aluminum  alloy  plate.  The  experimental 
data,  as  well  as  the  numerical  results  S  obtained  from  equation  (3)  are 
tabulated  in  Table  1.  From  Table  1,  it  can  be  concluded  that  the 
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numerical  values  S  calculated  using  equation  (3)  are  relatively  close 
to  the  average  values  SE  of  the  actual  measurement . 


Table  1  Comparison  between  the  calculated  S  and  the  tested  St 
in  a  2024-T3  aluminum-  alloy 


» 

I- 

i 

SKot 

MPa./m 

C  ft 

Rol 

C  8  3 

MPa«/m 

A  K,» 

C  9  3 

MPa«/  m 

3ia*rn.h*s%ia 

mm 

C  7  3  ( 3 

1  «•*«  0 

!  r.  _  AX’oj,-  A  Ktk 

|!  (l~RoQK, 

I.  I  | 

24.7 

0.13 

88.0 

3.4 

0.26 

1  0-2S 

I 

24.7 

0.21 

81.0 

2.8 

0.34 

0.32 

i  i 

24.7 

0.57 

38.0 

1.7 

0.67 

0.61 

Key:  1 — experimental  group  number;  2 — overload  cyclic  stress  ratio; 

3 — averagt.  measured  value  S£  of  the  ratio  of  slanted  cracks 
shear  type)  in  fractured  cracks;  4 — calculated  value 

III.  IMPROVEMENT  FOR  MAARSE  RETARDATION  MODEL 

Maarse  model  [10]  is  based  on  the  crack  closure  theory  and  its 
basic  equation  is 

where  Kjjp  is  the  opening  stress  intensity  factor  corresponding  to  the  crack 

opening  •  load  PQp.  Maarse  used  the  "residual  pressure  stress  crack 

negative  opening  displacement"  assumption  to  propose  an  engineering 

method  to  predict  P„  .  i.e.,  „ 

op 1  p  r _ tii  ct _ 

/  a .  \ 

(5) 

where  „  B(2*fV)l,t 

8  (  i  _v»)  (plane  strain) 

where  B  and  W  are  the  thickness  and  width  of  the  specimen  respectively; 

/ is  a  relevant  function  of  the  crack  length  and  specimen  geometric 

shape  in  the  expression  of  stress  intensity  factor;  (a  -a)  is  the 

s 

length  in  the  x  direction  in  the  residual  plastic  compression  zone; 

/fret  is  twice  the  measured  distance  from  the  crack  tip  to  the  plastic 
compression  boundary  along  the  Y  direction  (Figure  1). 

The  use  of  Maarse  model  to  estimate  the  effect  of  overload  retard¬ 
ation  has  two  inconvenient  points:  one  is  that  the  material  constants 
C*  and  n*  in  equation  (4)  must  first  be  treated  according  to  AKeff  and 
the  other  is  that  the  calculation  under  mixed  type  and  stress  conditions 

was  not  given  in  the  original  paper. 
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With  respect  to  the  above  two 
points,  we  attempted  to  make  sim¬ 
plifications  and  additions  for  the 
Maarse  model  in  this  paper. 


1.  Simplification  of  the  calcula¬ 
tion  of  C*  and  n* 


U  k  Vi r(amy 

«  w; 


( _I_*J  'll?  1  ( _  *  , 

■)  ZZi  ‘  "  1 


The  Paris  equation  to  calculate  V .  | 

the  rate  of  crack  growth  under  con-  - - - -  * 

stant  amplitude  is  - — - 

-j^=C(lK)m-c:Kw(  i  -/?»;••  (6) 

_  _  ,  ,  *"'*■  1  Residual  piastic  zone  geometry  and  location 

If  we  assume  that  the  crack  clo-  ... 

sure  stress  and  opening  stress  are  (  x  - ;  1 

0.0«8  ■)  ZZi  ;  =  1 

equal,  and  define  that  the  closure 
factor  C  is 

P  K  (7> 

/—p - K 

i  au  *  *  max 

then  the  Maarse  equation  (4)  used  to  calculate  da/dN  can  be  rewritten 
as : 

^-C*Cff«.(l-C,)r  (8) 

Based  on  the  experimental  results  in  [11],  equation  (6)  at  diff¬ 
erent  R  is  a  series  of  approximately  parallel  lines  with  difference 

dSL 

intercepts  in  the  jjj  AK  log-log  coordinate  system.  This  means  that 
the  change  of  the  exponent  n  is  very  small.  Hence,  it  can  be  assumed 
that 

( Q) 


In  fact,  for  the  same  condition  regardless  of  the  expression,  the 
crack  propagation  speed  is  the  same.  Then  equation  (6)  should  be  equal 
to  equation  (8):  l l  -C,)3* 

^_r.(  1  ~R  _ IjlB. _ T  (10) 


In  the  process  of  additive  calculation  of  crack  expansion,  it  is 


possible  to  let  the  initial  PQp  =  0.  The  PQp  in  other  stress  cycles 
can  be  computed  based  on  equation  (5). 

Thus,  after  the  simplification  through  equations  (9)  and  (10),  C* 
and  n*  can  be  calculated  using  the  material  constants  C  and  n  in  equa¬ 
tion  (6).  Consequently,  the  application  of  engineering  predictions  of 
crack  propagation  life  is  becoming  more  convenient. 

2.  The  extension  of  residual  plastic  zone 
compressible  ellipse  equation 

In  the  Maarse  model,  Rycg  an  important  parameter  in  the  calcula 
tion  of  PQp.  However,  in  the  original  paper  only  the  ellipse  equation 
determining  RycB  un<^er  plane  strain  conditions  was  given  which  limited 
the  application  of  this  model. 


In  this  paper  we  attempted  to  extend  the  original  ellipse  equation 
to  more  general  conditions,  including  plane  strain,  mixed  type  and  plane 
stress  conditions.  The  adopted  method  is  the  corrected  equation  of  the 
plastic  zone  dimension  described  in  the  previous  section  of  this  paper. 


i.e. , 


(11) 


a  is  calculated  based  on  equations  (3)  and  (2).  Thus,  the  ellipse  equa¬ 
tion  in  the  original  model  is  replaced  by  the  following  equation: 


1 


(12) 


IV.  CALCULATED  RESULTS  AND  DISCUSSIONS 


In  this  paper  the  overload  retardation  effect  was  calculated  using 
the  additive  method  in  combination  with  improved  Willenborg  model  and 
Maarse  model  to  predict  the  crack  propagation  behavior  in  T1-6A1-4V 
titanium  plate  overloading  equipment.  The  results  were  compared  with 
those  of  original  corresponding  models  and  experimental  results  [1], 
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The  relevant  data  for  the  computation  are  the  following: 

The  specimen  is  300x100x2  mm  of  the  center  crack  tensile  (CCT)  type, 
the  basic  cyclic  stress  ratio  1.6kN/l6.oicN=o.i» 

the  yield  stress  of  the  material  a,/-917.0MPa»  £.-115.4MPa»/mcl*\ ±K,k 
=6.5MPa»/m:‘S3  ;  the  material  constants  C  and  n  are  listed  in  Table  2 

(in  Paris  equation,  a2C>  13. 5 MPa^a):  Poisson  ratio  v-0.33. 


Table  2  Material  constants  in  FGG  for  Ti-6A1**4V 


ooircutational  equation 

c  ! 

n 

«  -cum-  ; 

1.341  xio*7 

2.515 

da  caio* 

dN~  (1  -  R)Kk- &K  i 

2.048 x  10'« 

2.311 

il 


Figure  2  shows  the  predicted 
retardation  cycle  number  and  mea¬ 
sured  number  Np  using  both  improved 
and  original  Willenborg  and  Maarse 

models  (overload  ratio  Qol-Po  .  IP _ ■»  t .  8) 

under  a  single  tensile  overload  ) . 
Figure  2  also  shows  the  variation 
of  a  before  and  after  the  improve¬ 
ment.  It  is  obvious  that  the  cor¬ 
rected  a  value  is  adjusted  automat¬ 
ically  (without  the  need  of  making 
a  selection  beforehand)  based  on  the 
variation  of  the  overload  level  KqL 
between  the  original  constant  values 
of  6^2.  This  advantage  may  be  more 
significant  in  the  prediction  of 
crack  propagation  life  under  complex 
spectrum  loads. 


Fi|.  2  Coefficient  «  and  number  of  delay  cycles  .V a 
versus  Kol  for  Tr-6Ai~4V  alloy  (Qj m  1.8) 


•1 


* 


It  is  also  possible  to  see  from  Figure  2  that  the  calculated  ND  cal 
values  of  plane  stress  (A  in  Figure  2)  and  plane  strain  (7  in  Figure  2) 
using  the  Willenborg  model  deviated  highly  from  the  actually  measured 
value  Nd  E.  The  former  was  shifted  toward  danger  and  the  latter  toward 
safety.  The  improved  Willenborg  model  predication  (in  Figure  2) 
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Table  3  Comparison  of  predicting  accuracy  of  modified  models 


with  original-  models 


•  •  ■£> 

Willenborg  US  (3> 

Willenborg  II 9 

Maarse  US 

<S> 

*  iSJS  69  Maarse  «S 

a’2® 

a-6fl 

X/Vdi  Cal 
3?o„ 

<»**is«i**« 

(P 

1.67 

0.53 

0.83 

0.37 

0.83 

Neal-  Sf 

(«*»*>  (§) 

V27.9SS 

-4.19K 

*  3.9966 

i 

i 

-  10.2SH  j  —  2.06% 

! 

i 


1 — error;  2 — Willenborg  model;  3 — plane  stress;  4 — plane  strain; 
5 — Improved  Willenborg  model;  6 — Maarse  model;  7 — Improved 
Maarse  model;  8 — (predicated  average  value  and  measured  average 
value);  9 — relative  error 


5  -li—ft  a  tad  Pft  vs  a  carves  is  *  single  tensile  overload  retardation  rone  for  Ti—SAl—  4V  alloy 
d.V 
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however,  was  located  in  between  the  above  results. 

The  calculated  a-N  curves  of  crack  propagation  under  multiple 
single  tensile  overload  based  on  the  improved  and  original  Willenborg 
and  Maarse  models  as  well  as  the  measured  a-N  curve  were  all  plotted 
in  Figure  3- 

Figure  4  also  shows  the  calculated  and  measured  values  of  the 
improved  and  original  Maarse  model  when  Qql  =  2.0. 

Summarizing  Figures  2,  3  and  4,  the  computational  accuracies  of 
the  above  models  before  and  after  the  improvement  are  listed  in  Table  3 

The  Ng  and  Ncal  in  Table  3  represent  the  measured  and  calculated 
number  of  load  cycles  necessary  to  propagate  a  crack  with  a  half  length 
a  *  11  mm  to  35  mm. 

Apparently,  from  Table  3,  the  accuracy  of  the  improved  Willenborg 
model  is  significantly  higher  than  that  of  the  original  model  and  its 
absolute  error  is  minimal.  Similarly,  the  improved  Maarse  model  has  a 
higher  prediction  accuracy  than  that  of  the  original  model  with  more 
satisfactory  results. 

For  further  comparison  and  exploration.  Figure  5  plots  the  retard¬ 
ation  characteristics  of  crack  propagation  under  single  tensile  over- 
da 

load:  -  a  vs.  PQp  -  a  curve.  From  Figure  5,  during  the  variation 

of  the  external  load  P  from  low-high-low,  the  variation  of  crack  open¬ 
ing  load  P _  as  described  by  the  Maarse  model  which  is  based  on  the 

op  da 

crack  closure  effect  and  the  corresponding  ^  is  shown.  When  the  peak 
load  P.~T  just  appears,  P_  immediately  drops  to  a  minimum  (caused  by 

vJii  Op 

the  split  tip  dullness  of  PQL)  and  then  it  gradually  rises  to  a  maximum 
(maximum  value  of  Ry0  corresponding  to  P0L)  and  begins  to  decline  until 
it  reaches  the  constant  amplitude  steady  PQp.  At  this  time  the  retard¬ 
ation  effect  disappears.  The  corresponding  variation  of  ^  is  opposite 

n  O 

to  that  of  PQp.  In  the  figure  it  also  indicates  that  the  curve 
calculated  based  on  the  improved  Maarse  model  approaches  the  experi¬ 
mental  curve  more  than  that  of  the  original  model.  Similarly,  the  retards- 


I 


retardation  model  based  on  the  closure  effect  can  better  describe  the 
five  stages  [1]  of  the  retardation  process  and  thus  are  conceptually 
better  from  a  physics  point  of  view. 

The  above  is  a  comparison  of  the  prediction  of  overload  retardation 
effect  for  the  material  T1-6A1-4V  under  two  overload  ratio  and  various 
crack  lengths.  As  for  the  application  of  other  materials  and  compli¬ 
cated  spectrum  loads,  further  verification  is  required. 

V.  CONCLUSIONS 

1.  The  plastic  zone  correction  equation  proposed  in  this  paper  is 


UJL)' 

«*\<w 


where 


s~(aa:-aa:,*)/(i  -r)k. 

The  corrected  plastic  zone  equation  can  more  closely  reflect  the 
characteristics  of  fatigue  crack  propagation  and  the  continuity  of 
stress  state  transformation.  Hence,  it  is  possible  to  be  suitable  for 
the  variable  stress  condition  under  variable  loads. 

2.  The  improved  Willenborg  model  and  Maarse  model  have  a  wider 
range  of  applicability  than  the  original  ones.  They  are  convenient  to 
use  and  the  accuracy  is  significantly  higher. 
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MODIFICATION  OF  THE  WILLENBORG  AND  MAARSE  MODELS 
AND  APPLICATION  TO  THE  LIFE  PREDICTION 
FOR  Ti— 6AI— 4V  ALLOY 

Zhang  Y cm gkui,  Gu  Mingda  and  Yan  Minggao 
(Institute  of  Aeronautical  Materials') 

Abstract 

The  prediction  accuracy  of  the  overload  retardation  effect  on  fatigue  crack 
growth  (FCG)  is  closely  related  to  the  formula  for  the  plastic  zone  adopted 
by  each  retardation  model.  Hence,  in  order  to  express  more  rationally  the 
behavior  of  fatxgue.-crack  growth  process  and  the  continuity  of  the  stress  state 
transformation  (from  plane  strain  to  plane  stress  mode),  the  coefficient  <*  in 
the  formula  for  the  plastic  zone  used  in  the  Willenborg  and  Maarse  models 
can  be  expressed  as  the  following] 

a=-rhs~ 

where 

A  K-tfC,h 

s — TT^rWT 

S  is  the  fraction  of  the  stress  part  in  a  fraccure  surface.  Meanwhile,  it  was 
found  that  the  material  constants  C*  and  a*  in  the  formula  of  the  Maarse 
model  can  be  substituted  by  C  and  a  obtained  from  the  constant  amplitude 
loading  tests. 

From  the  above  mentioned  modification  of  the  Willenborg  and  Maarse 
models,  it  is  recognized  that  the  modified  models  can  be  applied  more  success¬ 
fully  and  conveniently  to  the  life  prediction  for  TH3AI-4V  alloy  under  a 
single  tensile  or  a  series  of  single  tensile  overloads.  The  results  calculated  with 
the  modified  models  were  found  to  be  fairly  consistent  with  the  experimental 
data. 
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STUDY  ON  THE  MEASURING  TECHNIQUE  OF  THIN  BOUNDARY  LAYERS 

Nanjing  Aeronautical  Institute,  He  Zhongwei 

ABSTRACT 

This  paper  introduced  a  self-designed  and  fabricated  miniature 
driver  with  measuring  probes  for  measurement  of  thin  boundary  layers. 
The  driver  used  the  Chinese  made  Hwei  model  28BF001  stepmotor  as  the 
main  design  body.  It  transforms  the  angular  displacement  of  the  motor 
into  linear  displacement.  The  minimum  linear  velocity  is  0.00625  mm 
per  unit  pulse.  Its  working  travel  is  5.0  mm.  It  has  good  static  char' 
acteristlcs.  The  outside  diameter  is  28  mm  and  height  is  50  mm.  The 
structure  is  simple.  In  the  paper  two  types  of  thin  boundary  layer 
probes  were  introduced. 

The  results  of  the  use  of  it  in  a  supersonic  flow  showed  that  the 
driver  worked  reliably.  The  probe  did  not  show  any  vibration  and  was 
able  to  measure  thin  boundary  layer  velocity  distribution  with  accuracy 

In  the  present  conditions  of  the  dimensions  of  supersonic  tunnels 
in  our  country,  it  has  very  practical  meaning  [1]  to  design  a  miniature 
driver  with  boundary  layer  probe  to  measure  thin  boundary  layers.  This 
paper  introduces  a  self-designed  and  developed  miniature  driver  and 
probe. 

I.  MINIATURE  DRIVER  FOR  MEASUREMENT  OF  BOUNDARY  LAYERS 

The  main  body  of  the  driver  is  a  domestic  model  28BF001  step  motor 
The  electrical  step  angle  is  3°.  The  maximum  static  rotation  moment 
is  0.18  kg. cm.  The  outside  diameter  is  28  mm,  length  is  28  mm. 

The  design  principle  of  the  driver  is  simple.  The  structure  is 
compact.  A  precisely  machined  M  6.5  x  0.75  fine  screw-treaded  rod  made 
of  durable  copper  material  was  lightly  heat  pressed  on  the  electrical 
shaft.  The  treaded  rod  Is  installed  with  precisely  machined  screw-nuts 


When  the  treaded  rod  produced  an  angular  displacement  along  with  the 
electrical  shaft  due  to  an  external  electric  pulse  signal,  the  screw 
nut  in  the  very  smooth  square  groove  together  with  the  probe  installed 
on  its  top  was  constrained  to  move  up  or  down.  When  fabricating  the 
driver,  precision  machining  is  required.  The  threaded  rod  and  the  screw 
nut  must  be  matched  so  that  after  installation  their  friction  force 
and  spacing  are  very  small. 

Figure  1  represents  the  204-2  driver  structural  schematic  diagram. 
The  main  data  are  in  the  following:  Minimum  linear  velocity  is  0.00625 
mm  per  pulse,  working  travel  is  5.0  mm;  the  static  characteristics  of 
the  driver  are  shown  in  Figure  2.  Curves  B  and  C  (reproducibility)  are 
static  characteristic  curves  determined  by  the  powerful  toolmicroscope . 

A  is  the  theoretical  curve;  B  and  C  are  very  close  to  A. 

The  driver  is  controlled  directly  using  the  domestic  BQDI-001  driv¬ 
ing  power  source.  It  is  convenient  to  use.  It  is  possible  to  make 
corrections  through  the  use  of  that  driving  power  source,  the  movement 
velocity,  direction  and  position  of  the  probe  fixed  at  the  top  of  the 
driver  can  be  controlled  by  frequency  variation,  fixing  direction  and 
fixing  phase. 

Received  MHy  1981. 


Fij.  1  Schematic.. of  driver 
2©4-*2  for  mcaiurmf  the  boun¬ 
dary  layer 


II.  DESIGN  OF  THIN  BOUNDARY  LAYER  TOTAL 

PRESSURE  PROBE 

Through  the  analysis  of  the  possible  displace¬ 
ment  of  the  probe  in  supersonic  flow,  we  designed 
a  special  probe  (Figure  3)  from  actual  practical 
experience.  In  order  to  avoid  the  effect  of  the 
probe  lever  on  the  measured  total  pressure  valve 
let  us  take  L/d  >10,  and  1/L-- ~  .  Beginning 
with  strengthening  the  stiffness  of  the  probe,  the 
tip  of  the  probe,  is  made  into  an  arc  shape  and 
the  probe  level  is  made  into  a  loop  structure. 

The  front  extending  portion  of  the 
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Fi*.  2  Tie  itatie  chmr«et*ri*tie*  of  dm*r  204-2 

1 —  static  curve  determined  using  the  all  purpose  microscope; 

2 —  static  curve  determined  by  reproducibility;  3 — theoretical 
curve;  4 — electrical  step  angle  3°;  5 — spacing  of  screw  thread; 
6 — electric  pulse  counter;  7 — linear  displacement 


|  204-2-11  0.2  10.11 1.010.51 
204-2-3)0. 14l0.07|l.oio.5! 


S3  204-2*1  a 'S-'kH- fit 
3  Sizes  of  *h*  nnuQtiary  ^robei  2^1 

8 — solder;  9 — surface;  lo—  unit 


probe  is  7-8°  sharply  split.  The  inner  channel  of  the  probe  is  an 
expansion  type  to  facilitate  the  pressure  measurement  response  speed.  •  , 

In  order  to  ensure  that  the  lower  surface  of  the  probe  and  the 
surface  of  the  object  have  a  linear  contact,  the  lower  surface  of  the 
forward  extension  portion  of  the  probe  and  the  surface  of  the  object 
maintain  a  3°  inclination  angle.  The  outside  heights  of  the  two  types 
of  probes  are  0.2  mm  and  0.14  mm  and  the  inside  heights  are  0.1  mm  and 
0.07  mm.  The  inside  width  and  outside  width  at  the  1.0  mm  opening  of  the  probe 
are  0.8  mm  and  0.6  mm,  respectively.  The  top  of  the  probe  was  made  into 
shape  by  directly  beating  0.8  mm  a  stainless  steel  tube  until  red  hot. 

The  sizes  of  the  probe  opening  were  measured  using  a  standard  plug  cal-  ■•••] 

liper.  The  external  shape  of  the  two  probes  is  the  same.  The  external  ; 

shape  and  size  were  examined  optically.  The  sizes  of  the  opening  of  the  102~ 
two  probes  are  equivalent  to  those  of  the  boundary  layer  probes  used 
internationally  [2,3]. 


III.  APPLICATION  EXAMPLES 


In  the  single-support  point  semi-soft  wall  nozzle  supersonic  tunnel 
at  Nanjing  Aeronautical  Institute,  the  model  204-2  driver  with  boundary 
layer  probe  was  used  to  measure  the  boundary  layer  parameters  on  the 
second  cone  of  the  intake  of  a  certain  airplane.  Experimental  condi-  •  ' 

tions  were:  Incoming  flow  Mach  number  M„— 1.97,  2.10  ;  downward  inclin¬ 

ation  angle  of  the  nozzle  a -o'  .  The  pressure  measuring  system:  SYD-1 
transducer  and  XJ-100  circuit  measuring  meter  equipped  with  model  LS-5  ^ 

digital  recorder. 


When  the  boundary  layer  of  the  cone  was  measured,  the  driver  was 
installed  inside  the  cone  and  the  probe  reached  out  on  the  cone  surface 
(Figure  4).  During  the  measurement  of  boundary,  the  probe  moved  grad¬ 
ually  from  the  side  out  starting  from  the  surface  according  to  a  step 
of  10  pulses  per  stop.  The  match  between  the  probe  level  and  the  cone 
body  is  very  important.  Five  pressure  measuring  tubes  were  specially 
installed  inside  the  cone  chamber.  The  measured  results  reflected  that 
internal  pressure  variation  in  the  cone  chamber  was  very  slight  before 
and  during  the  experiment.  This  explains  that  the  boundary  layer  of 
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the  cone  surface  did  not  have  a  conical  inward  flow  (Figure  4).  In 
addition,  right  underneath  the  probe  opening,  a  0.6  mm  diameter  sta¬ 
tic  pressure  tube  was  installed  on  the  cone  surface.  It  is  electri¬ 
cally  isolated  from  the  cone.  In  the  experiment  the  static  pressure 
tube  carried  electricity.  When  the  probe  and  the  end  surface  of  the 
static  pressure  tube  come  into  contact,  an  indicator  light  comes  on. 
The  extinguished  light  was  a  symbol  for  the  starting  of  the  boundary 
layer. 


Fit.  4  Portion  of  driver  204-2  with  •  probe  ia  •  bicone. 


Key:  1 — electric  Insulation;  2— solder;  3 — to  transducer;  4 — T  press 
ure  tube  to  transducer;  5 — static  pressure  tube;  6 — size  80 
particle  size;  7 — 204-2  transducer 


Fit.  5  Schliem  photograph  of  tbock  wore 
oattcra  o«  the  hirrnie 


Experimental  observation  showed  that  the  driver  and  probe  could 
work  steadily  in  supersonic  flow.  Only  during  the  start-up  of  the 
tunnel,  because  of  the  instability  of  the  flow,  the  probe-surface  con¬ 
tact  indicator  light  twinkled  weakly.  However,  after  the  start-up  of 
the  tunnel,  the  indicator  light  stayed  on  constantly.  In  the  experi¬ 
ment  observation  was  made  especially  with  respect  to  the -motion  of  the 
probe.  On  the  frame  of  the  observation  window  of  the  tunnel,  a  0.25  mm 
diameter  steel  wire  was  fixed  in  the  transverse  direction  and  the  probe 
opening  was  placed  close  to  the  image  of  the  steel  wire.  Then,  the 
tunnel  was  started  and  after  repeated  observations,  the  probe  was  not 
found  to  move  relative  to  the  steel  wire  in  stable  supersonic  flows. 

No  probe  vibration  was  found  either. 

Schlieren  photograph  of  the  probe  and  shock  wave  pattern  on  the 
bicone  under  the  super  critical  condition  at  the  intake  at  M.-2.10,  a —o' 
is  shown  in  Figure  5- 

Figure  6  shows  the  boundary  layer  velocity  distribution  of  the  103 
second  conical  surface  of  the  center  cone  at  the  station  60  mm  from  the 
conical  tip.  The  boundary  layer  velocity  distribution  of  the  same 
cross-section  measured  by  204-2-1  and  204-2-2  probes  are  shown  in  the 
figure.  In  the  figure  y  represents  the  perpendicular  distance  from  the 
probe  center  to  the  surface;  6 — thickness  of  the  boundary  layer;  u — 
local  flow  speed  at  the  boundary  layer;  — boundary  layer  interface 

flow  velocity;  N — exponent. 

In  comparison  of  velocity  distribution,  the  probe  was  placed  as 
far  away  from  the  surface  as  possible  not  affecting  the  static  pressure 
underneath.  This  value  was  used  to  represent  the  static  pressure  of 
the  entire  cross-section.  As  for  speed  in  the  boundary  layer  above 
the  sonic  speed,  a  correction  must  be  made  based  on  positive  shock  wave 
loss  in  order  to  obtain  the  actual  speed  of  the  wave  front. 


u/U-*  <y/S>u* 

-  W-4.56 

-  *204-2-1  tSttM  1 

-  1204-2-21SW-3I. 

***»«*  <«4> 


Fig.  6  Boundary  layer  profiles  on  the  second  conical  surface  at  the  station  60®  ®  from  the  conical  tip 

Key:  1 — use  204-2-1  probe j  2 — use  2-4-2-2  probe;  3 — conical  tip 
has  turning  band  (Figure  4);  4 — 204-2-1  type;  5 — 204-2-2 
type;  6 — unit:  mm 

CONCLUSIONS 

Under  the  condition  M.-2.10,  1.9?  and  a  «o’  ,  for  a  1:10  scale 
Intake  second  conical  surface  of  a  certain  airplane,  experimental 
results  of  boundary  layer  verified: 

1.  Model  204-2  driver  can  completely  satisfy  the  requirements  of 
measuring  thin  boundary  layer.  The  minimum  linear  velocity  of  the 
driver  Is  0.00625  mm  per  pulse;  working  travel  is  5.0  mm.  It  has  good 
static  characteristics.  The  experimental  and  theoretical  values  are 
very  close.  It  works  reliably  and  its  structure  is  simple.  Its  success 
ful  miniaturization  has  provided  a  very  important  means  in  the  small 
scale  boundary  layer  control  study. 


2.  204-2  type  boundary  layer  probe  size  is  equivalent  to  that 

used  internally.  The  probe  has  excellent  rigidity.  It  shows  no  vibra¬ 
tion  in  high  speed  flow.  The  probe  can  work  steadily. 


W 


After  more  than  half  a  year  of  actual  use.  It  is  verified  that 
the  mentioned  driver  and  probe  for  measurement  of  boundary  layer  can 
completely  satisfy  the  requirements  in  measuring  thin  boundary  layers. 
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STUDY  ON  THE  MEASURING  TECHNIQUE  OF  THIN 
BOUNDARY  LAYERS 

He  Zhongivei 

( Nanjing  Aeronautical  Institute) 

Abstract 

A  type  of  miniature  driver  with  measuring  probe  204-2  designed  by  Our¬ 
selves  for  measurement  of  thin  boundary  layers  is  described.  In  the  design  the 
homemade  miniature  stepmoior  28BF001  serves  as  a  main  body.  The  driver  is 
simple  in  construction  (see  Fig.  l)  and  reliable  in  operation.  It  transforms  the 
angular  displacement  of  the  motor  into  linear  displacement.  Its  minimum  linear 
velocity  is  0.00625mm  per  unit  pulse,  and  its  working  travel  is  5.0mm.  The 
driver  with  an  outside  diameter  of  28mm  and  height  of  50mm  has  good  static 
characteristics  (Fig.  2).  It  is  found  that  the  experimental  curve  A,  C  extremely 
coincides  with  the  theoretical  curve  B. 

J  Two  types  of  measuring  probes  204-2  are  ftiso  presented  in  this  paper. 
The  outer  height  is  0.2mm  and  inner  height  0.1mm  for  one  probe  tip,  while 
for  another  they  are  0.14mm  and  0.07mm  respectively.  Both  tips  possess  the 
same  outer  width  of  1.0mm  and  the  inner  one  of  0.6mm.  The  configurations 
of  both  probes  are  the  same  (see  Fig.  3). 

The  driver  and  the  probes  presented  above  have  already  been  used  for 
measuring  the  boundary  layer  parameters  on  the  second  cone  of  a  supersonic 
intake  working  in  supercritical  conditions  (Figs  4  and  5).  The  experiments  were 
carried  out  in  our  supersonic  intake  tunnel  with  free-stream  Mach  numbers 
1.97  and  2.10.  The  boundary  layer  profiles  were  shown  in  Fig.  6. 

By  observing  the  movement  of  the  probes  during  the  experiments,  it  was 
proved  that rthe, driver  and  probes  worked  steadily  in  a  supersonic  flow  with¬ 
out  vibration*  and  the  thin  boundary  layer  profiles  measured  were  quite 
correct. 
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MICROMECHANICAL  ANALYSIS  OP  CREEP  CRACK  GROWTH 
ON  A  NICKEL  BASED  SUPERALLOY 


Beijing  Iron  and  Steel  Research  Institute 

Shen  Huiwang,  Gao  Zhentao,  Liu  Changfu  and  Cal  Qigung 

ABSTRACT 

The  stress  rupture  life  test  of  smooth  and  precracked  superalloy 
specimens  shows  that  a  small  pre-existing  crack  can  seriously  reduce 
the  stress  rupture  strength  and  life  of  superalloy  components. 

Scanning  electron  microscopy  and  metallograpbic  analyses  indicate 
that  the  crack  propagates  along  grain  boundaries  by  cavitation.  Based 
on  micromechanical  analysis  of  the  cavity  growth  and  coalescence  with 
the  main  crack  tip,  creep  crack  propagation  equations  were  derived  based 
on  creep  J  integral  parameter  J  and  crack  opening  displacement  (COD) 
rate  6.  Experimental  data  of  GH^A  superalloy  proved  the  above  fracture 
mechanical  analysis. 

I .  INTRODUCTION 

The  crack  formation  and  propagation  analyses  on  the  bottom  of  the 
outer  wheel  of  the  turbine  or  at  the  root  of  the  groove,  or  in  the 
blade  are  complicated  because  of  hot  corrosion  and  the  interaction  be¬ 
tween  strain  fatigue  and  creep.  The  high  temperature  zone  of  the  engine 
sometimes  has  (or  through  hot  corrosion  forms) small  cracks.  Under  the 
high  temperature  creep  conditions  these  small  cracks  interact  with  the 
grain  boundary  cavities  which  causes  these  pre-existing  small  cracks  to 
rapidly  propagate. 

In  the  past,  long  crack  specimens  were  used  in  creep  crack  growth 
experiments  which  were  limited  to  certain  creep  conditions.  However,  in 
high  temperature  components,  most  of  the  cracks  formed  by  welding  and 
hot  corrosion  have  much  smaller  dimensions  than  that  of  the  cross- 
section  of  the  components.  In  addition,  periodic  check-up  and  mainte¬ 
nance  also  would  not  allow  the  crack  to  grow  too  long.  From  the  view- 


point  of  design  life  and  fracture  control,  it  is  necessary  to  carry  out 
experimental  analysis  on  the  rupture  life  and  crack  propagation 
using  superalloy  specimens  with  pre-existing  small  cracks  under  total 
creep  yield  conditions.  This  paper  is  based  on  this  consideration. 

A  study  on  the  creep  crack  propagation  of  superalloy  was  performed  by 
understanding  the  mechanisms  of  grain  boundary  cavity  formation,  growth 
and  coalescence  due  to  creep  strain  during  each  flight  cycle  of  the 
engine  turbine. 

II.  EXPERIMENTAL  METHOD  AND  MATERIAL 

The  superalloy  (GH^A)  used  in  the  experiment  was  melted  in  a 
three-ton  electric  furnace  and  then  fabricated  into  material.  The 
major  chemical  components  of  the  alloy  are  Ni— 20^>Cc— 2.5<?&Ti— 0.79&AI— l.sq&Nb, 
The  heat  treatment  schedule  is  1080°Cx8  hours  solid  solution.  750°C 
xl6  hours.  The  constant  temperature  mechanical  characteristics  are: 
aB  -  1226,  a0>2  -  814,  S  -  31*,  *  -  35%. 

Received  January  1981. 

The  specimens  used  in  the  creep  test  were  cylindrical  and  plate 
specimens.  The  dimensions  of  cylindrical  specimens  are  <|>22xl00  (mm), 
and  4>10xl00  (mm)  and  the  plate  specimen  sizes  are  3x20x60  (mm)  and 
3x25x60  (mm). 

For  the  10  mm  diameter  cylindrical  specimen  and  plate  specimens, 
the  cracks  were  pre-fabricated  at  ambient  temperature  and  then  creep 
fracture  experiments  were  carried  out  at  high  temperature. 

The  fractured  surface  was  examined  using  the  British  made  SEM-S^-10 
for  observation. 

III.  ANALYSIS  OF  EXPERIMENTAL  RESULTS  AND  DISCUSSION 

From  the  stress  rupture  life  experiment  results  obtained  for 
smooth  and  pre-fabricated  cracked  specimens  of  superalloy  as  shown  in 
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Fi«.  1  Stre*  rapture  lif,  of  smooth  and  pm  racked  superaUoy 
specimen*  at  700*C 

Key:  1 — smooth  specimen;  2 — prefabricated  crack,  (mm);  3 — (MPa) 

4 — prefabricated  crack,  (mm);  5~ prefabricated  crack,  (mm); 
6 — life  (hrs) 


Figure  1,  it  is  seen  that  pre-existing  cracks  can  seriously  decrease  the 
stress  endurance  life  and  strength. 

Metallographic  analyses  of  the  fractured  surface  and  measured  sur¬ 
face  indicated  that  the  crack  propagated  along  the  grain  boundary 
(Figure  2).  Concentrated  spots  were  located  on  the  fractured  surface 
along  the  grain  boundary  (Figure  2(a)  which  showed  the  crack  propaga¬ 
tion  was  realized  through  the  growth  of  the  grain  boundary  cavities 
and  the  coalescence  with  the  main  crack.  Therefore,  it  is  necessary  to 
establish  the  law  of  creep  crack  propagation  along  the  grain  boundary 
based  on  the  microplastic  mechanical  analysis  of  grain  boundary  cavita¬ 
tion  growth  and  coalescence  with  the  main  crack  under  the  strain  field 
caused  by  the  creep  stress  in  front  of  the  crack.  In  this  paper,  the 
crack  opening  displacement  rate  and  creep  J  integral  parameter  J  are 
used  to  replace  stress  intensity  factor  K.  Because  under  high  tempera¬ 
ture  and  high  stress  conditions,  in  the  total  creep  yield  region  of 
existing  small  cracks,  the  creep  rate  is  high  which  quickly  relaxes  the 
elastic  stress  field,  hence  the  crack  tip  field  is  controlled  by  power 
multiplication  of  the  equation  of  creep  regularity.  For  material  with  a 
power  multiplication  creep  regularity  (e-  aaa)  ,  the  equations  of  mech¬ 
anics  and  the  solutions  are  completely  similar  to  those  of  pure  power 
solid  material  [2]  (where  t  is  creep  rate;  n  is  creep  exponent;  a  is 
a  material  constant  related  to  temperature;  a  is  stress).  For  a 


spherical  bole  with  an  instantaneous  radius  R  in  an  infinitely  large 
creeping  material,  its  radial  displacement  rate  R  can  be  expressed  as 

«*“*(««)  (1) 

where  £(«)  is  a  parameter  related  to  the  rigidity  of  the  material. 

For  a  pure  rigid  material,  equation  (1)  can  be  obtained  using  the 
similarity  analysis  method  according  to  the  literature  [2].  Let  us 
assume  that  the  rate  of  creep  cavitation  growth  in  the  crack  tip  region 
can  be  approximated  by  the  cavitation  creep  growth  rate  described  above 
It  is  worthwhile  noticing  that  only  when  the  cavity  is  very  small  rela¬ 
tive  to  the  space  where  the  cavity  exists  that  the  above  assumption  is 
correct.  Because  the  cavity  growth  life  is  mainly  spent  in  the  initial 
growth  period  when  the  cavity  is  still  small,  this  approximation  is  rea¬ 
sonable.  Then  we  get: 

where  a(x)  and  £(x)  can  be  any  arbitrary  stress  and  strain  rate  compo¬ 
nents  at  the  tip  of  the  crack. 


Under  steady  crack  propagation  conditions,  it  is  possible  to  cal¬ 
culate  the  crack  growth  rate  by  considering  that  the  crack  tip  moves 
toward  the  cavity  or  the. cavity  moves  toward  the  crack  tip.  If  A  is 
the  average  ditance  between  cavities,  under  the  condition  that  steady 

da  , '  in  the  process  that  the  cavity 


state  propagation  rate  is  b--^ 
approaches  the  crack  tip  at  x  « 


from  the  distance 


x  *  x  .  the  radius 
o 

of  the  cavity  grows  from  the  critical  radius  RQ  to  the  final  radius 
to  coalesce  with  the  crack  tip  which  ensures  the  steady  state 
propagation  of  cracks.  When  considering  the  cavity  growth  at  a  fixed 
point,  there  should  be  a  +  x  ■  constant.  Therefore, 

By  considering  dii(  x  x  at  )  •  -and  substituting  into  equation 


(2),  we  get 


and 


-ZLlL  r 


(3) 


i(*)dx 


•  lai7?r  J-f 

wbere  xQ  is  the  dimension  of  the  effect  stress  field  region  at  the 
crack  tip  which  is  the  inception  nucleus  of  the  cavity,  and  a  coordinate 
Increases  at  significant  rates  and  is  used  as  the  cut-off  distance  and 
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lower  limit  of  the  Integrals  of  the  right  hand  side  of  the  above  equa¬ 
tion.  Therefore,  for  materials  following  power  product  creep  law,  it  is 
possible  to  obtain  from  the  analytical  method  in  [2]  that: 

*(*>-(4-M‘V'  •)  <*) 

where  £(x)  is  the  arbitrary  strain  component  at  x.  Substitute  into 
equation  (3)  to  obtian  the  crack  propagation  law  described  by  the  crack 
opening  displacement  rate  parameter  6: 

■%-Bh  (5) 

20 

can  be  experimentally  determined  by  creep  fracture  experiments. 


Using  the  singularity  relation 
k(  x ) *  a o ,  (where 

I(0,  b)  )  is  the  angular  factor  para¬ 
meter  ofe(fi,  B)e-o),  the  crack  propa¬ 
gation  law  described  by  the  creep  J 
integral  parameter  J  can  also  be  ob¬ 


tained  by  substituting  it  into  equa¬ 


tion  (3). 


Cj«*» 


(6) 


where  c-(, (.)/ln^)j^,(_2_)-^ = o.  . 

This  relation  agrees  with  a  large 

amount  of  experimental  data  provided  in  [6] 


Fig.  3  Growth  oi  cavities  formed  alone 
grain  boundaries  under  crack  tip  creep 
stress  strain  field  and  their  coalescence 
with  the  main  crack  tip 


Sbib  and  Hutchinson  [7]  gave  the  finite  element  results  of  power 
product  law  creep  material  plane  stress  crack  tensile  plate  5: 

6-«K.(-2-.  «)«(-^rr)  (7) 

where  b  is  the  width  of  the  plate.  Substituting  equation  (7)  into  (5), 

we  get  da  m  n/  bo 

dl  \  b  —  a 
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v  0  '  which  can  be  treated  as  a  constant  in  engineering 

applications.  Equation  (8)  begins  to  integrate  from  the  initial  crack 
depth  a^  to  the  final  depth  a^.  to  obtain 


(9) 


At  700°C,  it  has  been  determined  that  the  creep  exponent  n  =  10  and 


,-26 


creep  crack  propagation  coefficient  D  =  0.70x10  for  superalloy  crack 
tensile  plate.  The  corresponding  experimental  results  are  shown  in 
Figure  1.  The  two  solid  lines  on  the  very  bottom  were  curves  obtained 


from  the  numerical  inttgration  of  equation  09)  by  using  n  =  10,  D  = 


0.70x10 
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*  2  mm,  a^-«  k  mm  and  af  *  10  mm.  The  points  in  the 


figure  are  the  actual  points  of  measurement  which  completely  coincide 
with  the  calculated  values.  Although  the  above  are  experimental  results 
obtained  from  thin  plates  under  creep  conditions,  it  is  possible  to 
extend  to  the  crack  growth  law  of  small  cracks  under  total  creep  yield 


thin  plate3.  When  — •><)  ,  the  integral  gives  the  result 


a  *«<expC0.70xi(rMo1,of  *  However,  for  the  small  surface  cracks  frequently 
encountered  in  engine  components,  these  cracks  are  not  located  in  the 
plane  stress  and  durable  band  creep  regions  but  primarily  in  the  plane 
strain  and  total  creep  region.  Hence,  a  stress  rupture  experiment 
(Figure  1)  was  performed  on  a  cylindrical  tensile  sample  with  prefab¬ 
ricated  surface  crack  at  a  depth  of  0.5-0. 8  mm  under  plane  strain  con¬ 
ditions.  Under  tough  zone  creep  conditions,  due  to  the  complexity  of 
geometry  involved  with  the  cylindrical  specimen  under  force,  the  anal¬ 
ysis  of  life  and  propagation  rate  of  the  surface  crack  becomes  rela¬ 
tively  more  difficult.  Therefore,  it  is  necessary  to  limit  the  experi¬ 
ment  under  total  creep  conditions.  It  is  to  say  that  an  interrupted 
stress  endurance  test  was  performed  using  tensile  rod  with  pre-existing 
cracks  at  0.5-0. 8  mm  depth  to  make  the  final  crack  depth  af  not  exceed 
2.5  mm.  Therefore,  under  total  creep  conditions  6 -2nQa(aom) 

obtained  from  [2]  can  be  plugged  into  equation  (5)  and  after  integration 


we  get 


Do"*-  ln- 


Oy 


(10) 


where  D-2-tQaB  j  q — crack  shape  factor 


3 


At  700°C,  under  plane  strain  conditions,  the  creep  exponent  of 
the  superalloy  was  experimentally  determined  to  be  n  *  10.  Under  total 
creep  yield  conditions,  the  same  result  as  equation  (10)  was  obtained 
by  substituting  j=2*Q*aJ odi™  into  equation  (6)  where  D- 

c[2n  Q*a~  *_  -j«~+T  .  The  photograph  in  Figure  2(b)  is  a  sideview  of  a 

metallographic  crack  from  the  interrupted  stress  experiment.  The  cut¬ 
ting  depth  is  0.72  mm.  The  fatigue  crack  depth  across  the  grain  is 
0.1  mm  measured  from  the  picture.  The  initial  crack  depth  a^  =  0.82 
mm.  The  creep  crack  along  the  grain  boundary  in  the  photograph  is  0.4 
mm.  Therefore,  a^  =  1.22  mm.  The  experimental  stress  a  =  373  MPa  and 
period  t  *  153  hours.  Putting  the  data  into  the  equations,  the  creep 

crack  propagation  coefficient  under  plane  strain  conditions  was  calcu- 

1  72 

lated  to  be  D~\ti-~r/Z7zl'x.  153-0.50 xio*Jt  Consequently,  the  crack  propaga¬ 
tion  formula,  a  —aiexpCO.SOxKr1*^*/)  is  obtained  under  plane  strain 

and  total  creep  conditions.  It  is  worthwhile  noticing  that  the  crack 
propagation  rate  in  this  paper  is  two  orders  of  magnitude  lower  under 
plane  strain  conditions  than  that  of  plane  stress  conditions.  SEM  anal¬ 
ysis  of  fractured  surface  showed  that  the  cavity  size  and  spacing  in 
the  middle  region  under  plane  strain  are  much  smaller  than  those  of 
plane  stress  specimens  (Figure  4).  This  explains  that  the  cavity  pro¬ 
pagation  rate  and  creep  rate  at  the  crack  tip  are  lower  under  plane 
strain  conditions.  Therefore,  before  the  coalescence  of  cavities  more 
cavitation  nuclei  were  formed.  Experimental  results  showed  that  under 
plane  strain  and  total  creep  conditions  small  surface  cracks  are  not 
the  same  as  those  penetrating  cracks  of  thin  plates  under  plane  stress 
and  durable  band  creep  conditions.  Hence,  in  order  to  more  precisely 
predict  the  life  and  crack  propagation  rate  of  engine  components  with 
pre-existing  surface  cracks,  it  is  necessary  to  carry  out  stress  endur¬ 
ance  life  and  crack  propagation  rate  analysis  and  test  using  specimens 
with  pre-existing  small  surface  cracks  under  plane  strain  and  total 
creep  conditions. 

IV  CONCLUSIONS 

Pre-existing  small  cracks  in  superalloy  components  can  seriously 
decrease  stress  rupture  life  and  strength. 
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Fig.  4  Cavities  formed  along  grain  boundaries  it  700*C  SEM  micrographs  of  fracture  surfaces  x  500 

<a)  Plane  straini  (5)  Plane  stress. 

On  the  basis  of  micromecbanical  analysis,  using  crack  opening 
displacement  rate  and  creep  J  integral  to  evaluate,  fracture  mechan¬ 
ical  analysis  has  been  performed  for  power  product  creep  material  crack 
tensile  specimens.  It  was  discovered  that  the  creep  propagation  rate 
is  two  orders  of  magnitude  lower  in  plane  strain  conditions  than  that 
in  plane  stress  conditions. 
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MICROMECHANICAL  ANALYSIS  OF  CREEP  CRACK 
GROWTH  ON  A  NICKEL-BASE  SUPERALLOY 

Shen  Huiwang,  Gao  Zhentao,  Liu  Changfu,  Cai  Qigung 
(.Central  Iron  and  Sieel  Research  Inslilute) 

Abstract 

Experiment  on  stress  rupture  life  of  smooth  and  precracked  superalloy 
specimens  shows  that  small  preexisting  crack  seriously  reduces  the  stress  rup¬ 
ture  strength  and  life  of  superalloy  components. 

SEM  and  metallographic  analysis  indicate  that  the  crack  grows  along  grain 
boundaries  by  cavitation.  Based  on  micromechanical  analysis  of  the  cavity 
growth  and  coalescence  with  the  main  crack  tip*  creep  crack  growth  formulas 
were  derived  by  using  creep  /  integral  parameter  j  and  COD  rate  parameter  6. 
Experimental  data  on  GH»A  superalloy  have  verified  the  above  fracture  mecha¬ 
nical  analysis. 


THE  "VIBRATION  THEORY  AND  APPLIED  SCIENCE  CONFERENCE" 


Jointly  held  by  China  Aeronautical  and  Astronautical  Society 
and  China  Mechanics  Society 


The  China  Aeronautical  and  Astronautical  Society  and  the  China 
Mechanics  Society  held  the  "Vibration  Theory  and  Applied  Science  Con¬ 
ference"  on  December  10-16,  1981  in  Quinming.  The  formal  representatives 
and  the  attending  representatives  were  168  in  total.  They  came  from  20 
provinces  and  cities  and  8l  organizations.  There  were  115  papers  pre¬ 
sented.  The  meeting  reviewed  the  research  results  of  vibration  theory 
and  its  application  in  our  country. 

Some  of  the  experts  attending  the  meeting  gave  technical  reports. 

As  examples,  papers  such  as  "Coupled  mechanical  problems  of  fluid  and 
structures"  by  Professor  Tu  Ching  hua,  "Identification  of  time-space 
parameters  and  on-board  reduction  methods"  by  Professor  Huang  Wenhu, 
"Application  of  vibration  theory  in  mechanical  fields"  by  Professor 
Chu  Weitu  and  "Structural  mechanics  of  space  craft"  by  Professor  Hu 
Haichan,  etc.,  benefited  the  audience. 


The  conference  was  divided  into  four  topics  of  exchange  informa¬ 
tion;  viz.,  simulation  combination  and  characteristics,  parametric 
identification  and  measurement  technique,  aerodynamic  elastic  on-board 
vibration  and  its  effect;  and  rotor  dynamics,  vibration  reduction  and 
isolation,  etc.  After  the  exchange,  all  the  representatives  believed 
that  vibration  theory  has  accomplished  some  results  in  the  application 
to  engineering  problems  in  our  country  during  recent  years.  By  absorb¬ 
ing  new  theories,  new  technologies  and  new  results,  this  science  has 
been  developed  and  experimental  techniques  and  measuring  equipment  have 
been  improved  significantly.  In  the  conference  the  22nd  Meeting  of 
Structural  Mechanics  and  Materials  in  the  US  was  introduced.  Further¬ 
more,  the  9th  International  and  1st  and  2nd  National  Non-linear 
Vibration  Meetings  were  also  introduced.  Hence,  the  status  of  vibra¬ 
tional  study  in  and  out  of  the  country  was  further  discussed. 
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After  exchanging  ideas,  the  direction  of  future  research  was 
explored.  It  was  proposed  that  the  trend  of  development  should  be  care 
fully  noticed.  The  studies  on  flow-solid  interaction,  active  control 
technology  and  time-space  method  in  conjunction  with  theoretical 
study,  engineering  applications  and  experimental  techniques  should  be 
further  strengthened  concurrently. 

During  the  conference,  the  China  Aeronautical  and  Astronautical 
Society  and  the  China  Mechanical  Society  held  meetings  to  discuss  the 
formation  of  vibration  special  groups  and  the  content  of  future  aca¬ 
demic  activities,  respectively. 

The  next  national  vibration  theory  will  also  be  held  jointly  by 
the  two  societies,  tentatively  in  1984. 
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